diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h
index d0474c8344..ee3419fa22 100644
--- a/include/openmc/hdf5_interface.h
+++ b/include/openmc/hdf5_interface.h
@@ -345,7 +345,7 @@ read_dataset(hid_t obj_id, const char* name, Position& r, bool indep=false)
 }
 
 template <typename T, std::size_t N>
-void read_dataset_as_shape(hid_t obj_id, const char* name,
+inline void read_dataset_as_shape(hid_t obj_id, const char* name,
                            xt::xtensor<T, N>& arr, bool indep=false)
 {
   hid_t dset = open_dataset(obj_id, name);
@@ -367,7 +367,7 @@ void read_dataset_as_shape(hid_t obj_id, const char* name,
 
 
 template <typename T, std::size_t N>
-void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
+inline void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
      bool must_have=false)
 {
   if (object_exists(obj_id, name)) {
diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h
index 1b4062da93..313677f6cf 100644
--- a/include/openmc/mgxs.h
+++ b/include/openmc/mgxs.h
@@ -110,7 +110,10 @@ class Mgxs {
     //!
     //! @param xs_id HDF5 group id for the cross section data.
     //! @param temperature Temperatures to read.
-    Mgxs(hid_t xs_id, const std::vector<double>& temperature);
+    //! @param num_group number of energy groups
+    //! @param num_delay number of delayed groups
+    Mgxs(hid_t xs_id, const std::vector<double>& temperature,
+        int num_group, int num_delay);
 
     //! \brief Constructor that initializes and populates all data to build a
     //!   macroscopic cross section from microscopic cross section.
@@ -119,8 +122,11 @@ class Mgxs {
     //! @param mat_kTs temperatures (in units of eV) that data is needed.
     //! @param micros Microscopic objects to combine.
     //! @param atom_densities Atom densities of those microscopic quantities.
+    //! @param num_group number of energy groups
+    //! @param num_delay number of delayed groups
     Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
-         const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities);
+         const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
+         int num_group, int num_delay);
 
     //! \brief Provides a cross section value given certain parameters
     //!
diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h
index bf72d675bd..afd183e87e 100644
--- a/include/openmc/mgxs_interface.h
+++ b/include/openmc/mgxs_interface.h
@@ -12,36 +12,71 @@
 namespace openmc {
 
 //==============================================================================
-// Global variables
+// Global MGXS data container structure
 //==============================================================================
 
+class MgxsInterface {
+public:
+
+  MgxsInterface() = default;
+
+  // Construct from path to cross sections file, as well as a list
+  // of XS to read and the corresponding temperatures for each XS
+  MgxsInterface(const std::string& path_cross_sections,
+                const std::vector<std::string> xs_to_read,
+                const std::vector<std::vector<double>> xs_temps);
+
+  // Does things to construct after the nuclides and temperatures to
+  // read have been specified.
+  void init();
+
+  // Set which nuclides and temperatures are to be read
+  void set_nuclides_and_temperatures(std::vector<std::string> xs_to_read,
+                                     std::vector<std::vector<double>> xs_temps);
+
+  // Add an Mgxs object to be managed
+  void add_mgxs(hid_t file_id, const std::string& name,
+         const std::vector<double>& temperature);
+
+  // Reads just the header of the cross sections file, to find
+  // min & max energies as well as the available XS
+  void read_header(const std::string& path_cross_sections);
+
+  // Calculate microscopic cross sections from nuclide macro XS
+  void create_macro_xs();
+
+  // Get the kT values which are used in the OpenMC model
+  std::vector<std::vector<double>> get_mat_kTs();
+
+  int num_energy_groups_;
+  int num_delayed_groups_;
+  std::vector<std::string> xs_names_; // available names in HDF5 file
+  std::vector<std::string> xs_to_read_; // XS which appear in materials
+  std::vector<std::vector<double>> xs_temps_to_read_; // temperatures used
+  std::string cross_sections_path_; // path to MGXS h5 file
+  std::vector<Mgxs> nuclides_;
+  std::vector<Mgxs> macro_xs_;
+  std::vector<double> energy_bins_;
+  std::vector<double> energy_bin_avg_;
+  std::vector<double> rev_energy_bins_;
+  std::vector<std::vector<double>> nuc_temps_; // all available temperatures
+};
+
 namespace data {
+  extern MgxsInterface mg;
+}
 
-extern std::vector<Mgxs> nuclides_MG;
-extern std::vector<Mgxs> macro_xs;
-extern int num_energy_groups;
-extern int num_delayed_groups;
-extern std::vector<double> energy_bins;
-extern std::vector<double> energy_bin_avg;
-extern std::vector<double> rev_energy_bins;
+// Puts available XS in MGXS file to globals so that when
+// materials are read, the MGXS specified in a material can
+// be ensured to be present in the available data.
+void put_mgxs_header_data_to_globals();
 
-} // namespace data
+// Set which nuclides and temperatures are to be read on
+// mg through global data
+void set_mg_interface_nuclides_and_temps();
 
-//==============================================================================
-// Mgxs data loading interface methods
-//==============================================================================
-
-void read_mgxs();
-
-void
-add_mgxs(hid_t file_id, const std::string& name,
-     const std::vector<double>& temperature);
-
-void create_macro_xs();
-
-std::vector<std::vector<double>> get_mat_kTs();
-
-void read_mg_cross_sections_header();
+// After macro XS have been read, materials can be marked as fissionable
+void mark_fissionable_mgxs_materials();
 
 //==============================================================================
 // Mgxs tracking/transport/tallying interface methods
diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h
index a8aed5d17c..3e9fc51cca 100644
--- a/include/openmc/xsdata.h
+++ b/include/openmc/xsdata.h
@@ -22,6 +22,7 @@ namespace openmc {
 class XsData {
 
   private:
+
     //! \brief Reads scattering data from the HDF5 file
     void
     scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
@@ -61,6 +62,9 @@ class XsData {
     void
     fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang);
 
+    //! Number of energy and delayed neutron groups
+    size_t n_g_, n_dg_;
+
   public:
 
     // The following quantities have the following dimensions:
@@ -98,7 +102,10 @@ class XsData {
     //! @param scatter_format The scattering representation of the file.
     //! @param n_pol Number of polar angles.
     //! @param n_azi Number of azimuthal angles.
-    XsData(bool fissionable, int scatter_format, int n_pol, int n_azi);
+    //! @param n_groups Number of energy groups.
+    //! @param n_d_groups Number of delayed neutron groups.
+    XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
+           size_t n_groups, size_t n_d_groups);
 
     //! \brief Loads the XsData object from the HDF5 file
     //!
diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp
index 3c29dd8e6d..6245bca9f9 100644
--- a/src/cross_sections.cpp
+++ b/src/cross_sections.cpp
@@ -155,7 +155,8 @@ void read_cross_sections_xml()
   if (settings::run_CE) {
     read_ce_cross_sections_xml();
   } else {
-    read_mg_cross_sections_header();
+    data::mg.read_header(settings::path_cross_sections);
+    put_mgxs_header_data_to_globals();
   }
 
   // Establish mapping between (type, material) and index in libraries
diff --git a/src/initialize.cpp b/src/initialize.cpp
index fd524a2de9..460db24363 100644
--- a/src/initialize.cpp
+++ b/src/initialize.cpp
@@ -260,8 +260,9 @@ void read_input_xml()
       read_ce_cross_sections(nuc_temps, thermal_temps);
     } else {
       // Create material macroscopic data for MGXS
-      read_mgxs();
-      create_macro_xs();
+      set_mg_interface_nuclides_and_temps();
+      data::mg.init();
+      mark_fissionable_mgxs_materials();
     }
     simulation::time_read_xs.stop();
   }
diff --git a/src/material.cpp b/src/material.cpp
index 9a1e206ac8..a5bc3dc138 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -368,7 +368,7 @@ void Material::normalize_density()
     // determine atomic weight ratio
     int i_nuc = nuclide_[i];
     double awr = settings::run_CE ?
-      data::nuclides[i_nuc]->awr_ : data::nuclides_MG[i_nuc].awr;
+      data::nuclides[i_nuc]->awr_ : data::mg.nuclides_[i_nuc].awr;
 
     // if given weight percent, convert all values so that they are divided
     // by awr. thus, when a sum is done over the values, it's actually
@@ -388,7 +388,7 @@ void Material::normalize_density()
     for (int i = 0; i < nuclide_.size(); ++i) {
       int i_nuc = nuclide_[i];
       double awr = settings::run_CE ?
-        data::nuclides[i_nuc]->awr_ : data::nuclides_MG[i_nuc].awr;
+        data::nuclides[i_nuc]->awr_ : data::mg.nuclides_[i_nuc].awr;
       sum_percent += atom_density_(i)*awr;
     }
     sum_percent = 1.0 / sum_percent;
@@ -726,7 +726,7 @@ void Material::init_bremsstrahlung()
 void Material::init_nuclide_index()
 {
   int n = settings::run_CE ?
-    data::nuclides.size() : data::nuclides_MG.size();
+    data::nuclides.size() : data::mg.nuclides_.size();
   mat_nuclide_index_.resize(n);
   std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), C_NONE);
   for (int i = 0; i < nuclide_.size(); ++i) {
@@ -994,11 +994,11 @@ void Material::to_hdf5(hid_t group) const
   } else {
     for (int i = 0; i < nuclide_.size(); ++i) {
       int i_nuc = nuclide_[i];
-      if (data::nuclides_MG[i_nuc].awr != MACROSCOPIC_AWR) {
-        nuc_names.push_back(data::nuclides_MG[i_nuc].name);
+      if (data::mg.nuclides_[i_nuc].awr != MACROSCOPIC_AWR) {
+        nuc_names.push_back(data::mg.nuclides_[i_nuc].name);
         nuc_densities.push_back(atom_density_(i));
       } else {
-        macro_names.push_back(data::nuclides_MG[i_nuc].name);
+        macro_names.push_back(data::mg.nuclides_[i_nuc].name);
       }
     }
   }
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index 1fca663080..84381fe77d 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -24,18 +24,6 @@
 
 namespace openmc {
 
-//==============================================================================
-// Global variables
-//==============================================================================
-
-namespace data {
-
-// Storage for the MGXS data
-std::vector<Mgxs> nuclides_MG;
-std::vector<Mgxs> macro_xs;
-
-} // namespace data
-
 //==============================================================================
 // Mgxs base-class methods
 //==============================================================================
@@ -53,8 +41,6 @@ Mgxs::init(const std::string& in_name, double in_awr,
   kTs = xt::adapt(in_kTs);
   fissionable = in_fissionable;
   scatter_format = in_scatter_format;
-  num_groups = data::num_energy_groups;
-  num_delayed_groups = data::num_delayed_groups;
   xs.resize(in_kTs.size());
   is_isotropic = in_is_isotropic;
   n_pol = in_polar.size();
@@ -284,7 +270,10 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
 
 //==============================================================================
 
-Mgxs::Mgxs(hid_t xs_id, const std::vector<double>& temperature)
+Mgxs::Mgxs(hid_t xs_id, const std::vector<double>& temperature,
+    int num_group, int num_delay) :
+  num_groups(num_group),
+  num_delayed_groups(num_delay)
 {
   // Call generic data gathering routine (will populate the metadata)
   int order_data;
@@ -299,7 +288,8 @@ Mgxs::Mgxs(hid_t xs_id, const std::vector<double>& temperature)
 
   // Load the more specific XsData information
   for (int t = 0; t < temps_to_read.size(); t++) {
-    xs[t] = XsData(fissionable, final_scatter_format, n_pol, n_azi);
+    xs[t] = XsData(fissionable, final_scatter_format, n_pol, n_azi,
+                   num_groups, num_delayed_groups);
     // Get the temperature as a string and then open the HDF5 group
     std::string temp_str = std::to_string(temps_to_read[t]) + "K";
     hid_t xsdata_grp = open_group(xs_id, temp_str.c_str());
@@ -317,7 +307,10 @@ Mgxs::Mgxs(hid_t xs_id, const std::vector<double>& temperature)
 //==============================================================================
 
 Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
-     const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities)
+     const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
+     int num_group, int num_delay) :
+  num_groups(num_group),
+  num_delayed_groups(num_delay)
 {
   // Get the minimum data needed to initialize:
   // Dont need awr, but lets just initialize it anyways
@@ -340,7 +333,7 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
   // Create the xs data for each temperature
   for (int t = 0; t < mat_kTs.size(); t++) {
     xs[t] = XsData(in_fissionable, in_scatter_format, in_polar.size(),
-      in_azimuthal.size());
+      in_azimuthal.size(), num_groups, num_delayed_groups);
 
     // Find the right temperature index to use
     double temp_desired = mat_kTs[t];
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index 29e5e8b363..1b0f2675a2 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -18,42 +18,53 @@
 
 namespace openmc {
 
-//==============================================================================
-// Global variable definitions
-//==============================================================================
-
-namespace data {
-
-int num_energy_groups;
-int num_delayed_groups;
-std::vector<double> energy_bins;
-std::vector<double> energy_bin_avg;
-std::vector<double> rev_energy_bins;
-
-} // namesapce data
-
 //==============================================================================
 // Mgxs data loading interface methods
 //==============================================================================
 
-void read_mgxs()
+namespace data {
+  MgxsInterface mg;
+}
+
+MgxsInterface::MgxsInterface(const std::string& path_cross_sections,
+                             const std::vector<std::string> xs_to_read,
+                             const std::vector<std::vector<double>> xs_temps)
 {
+  read_header(path_cross_sections);
+  set_nuclides_and_temperatures(xs_to_read, xs_temps);
+  init();
+}
+
+void MgxsInterface::set_nuclides_and_temperatures(
+    std::vector<std::string> xs_to_read,
+    std::vector<std::vector<double>> xs_temps)
+{ 
+  // Check to remove all duplicates
+  xs_to_read_ = xs_to_read;
+  xs_temps_to_read_ = xs_temps;
+  if (xs_to_read_.size() != xs_temps.size())
+    fatal_error("The list of macro XS temperatures to read does not "
+                "correspond in length to the number of XS names. ");
+}
+
+void MgxsInterface::init()
+{
+
+  // Check that at least some data was set to be read
+  if (xs_to_read_.size() == 0)
+    warning("No MGXS nuclides were set to be read.");
+
   // Check if MGXS Library exists
-  if (!file_exists(settings::path_cross_sections)) {
+  if (!file_exists(cross_sections_path_)) {
     // Could not find MGXS Library file
-    fatal_error("Cross sections HDF5 file '" + settings::path_cross_sections +
+    fatal_error("Cross sections HDF5 file '" + cross_sections_path_ +
       "' does not exist.");
   }
 
   write_message("Loading cross section data...", 5);
 
-  // Get temperatures
-  std::vector<std::vector<double>> nuc_temps(data::nuclide_map.size());
-  std::vector<std::vector<double>> dummy;
-  get_temperatures(nuc_temps, dummy);
-
   // Open file for reading
-  hid_t file_id = file_open(settings::path_cross_sections, 'r');
+  hid_t file_id = file_open(cross_sections_path_, 'r');
 
   // Read filetype
   std::string type;
@@ -73,38 +84,18 @@ void read_mgxs()
 
   // ==========================================================================
   // READ ALL MGXS CROSS SECTION TABLES
-
-  std::unordered_set<std::string> already_read;
-
-  // Build vector of nuclide names
-  std::vector<std::string> nuclide_names(data::nuclide_map.size());
-  for (const auto& kv : data::nuclide_map) {
-    nuclide_names[kv.second] = kv.first;
-  }
-
-  // Loop over all files
-  for (const auto& mat : model::materials) {
-    for (int i_nuc : mat->nuclide_) {
-      std::string& name = nuclide_names[i_nuc];
-
-      if (already_read.find(name) == already_read.end()) {
-        add_mgxs(file_id, name, nuc_temps[i_nuc]);
-        already_read.insert(name);
-      }
-
-      if (data::nuclides_MG[i_nuc].fissionable) {
-        mat->fissionable_ = true;
-      }
-    }
-  }
+  for (unsigned i_nuc=0; i_nuc<xs_to_read_.size(); ++i_nuc)
+    add_mgxs(file_id, xs_to_read_[i_nuc], xs_temps_to_read_[i_nuc]);
 
   file_close(file_id);
+
+  create_macro_xs();
 }
 
 //==============================================================================
 
 void
-add_mgxs(hid_t file_id, const std::string& name,
+MgxsInterface::add_mgxs(hid_t file_id, const std::string& name,
   const std::vector<double>& temperature)
 {
   write_message("Loading " + std::string(name) + " data...", 6);
@@ -118,13 +109,14 @@ add_mgxs(hid_t file_id, const std::string& name,
                 + "provided MGXS Library");
   }
 
-  data::nuclides_MG.emplace_back(xs_grp, temperature);
+  nuclides_.emplace_back(xs_grp, temperature, num_energy_groups_,
+      num_delayed_groups_);
   close_group(xs_grp);
 }
 
 //==============================================================================
 
-void create_macro_xs()
+void MgxsInterface::create_macro_xs()
 {
   // Get temperatures to read for each material
   auto kTs = get_mat_kTs();
@@ -140,24 +132,25 @@ void create_macro_xs()
       std::vector<double> atom_densities(mat->atom_density_.begin(),
         mat->atom_density_.end());
 
-      // Build array of pointers to nuclides_MG's Mgxs objects needed for this
+      // Build array of pointers to nuclides's Mgxs objects needed for this
       // material
       std::vector<Mgxs*> mgxs_ptr;
       for (int i_nuclide : mat->nuclide_) {
-        mgxs_ptr.push_back(&data::nuclides_MG[i_nuclide]);
+        mgxs_ptr.push_back(&nuclides_[i_nuclide]);
       }
 
-      data::macro_xs.emplace_back(mat->name_, kTs[i], mgxs_ptr, atom_densities);
+      macro_xs_.emplace_back(mat->name_, kTs[i], mgxs_ptr, atom_densities,
+          num_energy_groups_, num_delayed_groups_);
     } else {
       // Preserve the ordering of materials by including a blank entry
-      data::macro_xs.emplace_back();
+      macro_xs_.emplace_back();
     }
   }
 }
 
 //==============================================================================
 
-std::vector<std::vector<double>> get_mat_kTs()
+std::vector<std::vector<double>> MgxsInterface::get_mat_kTs()
 {
   std::vector<std::vector<double>> kTs(model::materials.size());
 
@@ -186,61 +179,112 @@ std::vector<std::vector<double>> get_mat_kTs()
 
 //==============================================================================
 
-void read_mg_cross_sections_header()
+void MgxsInterface::read_header(const std::string& path_cross_sections)
 {
+  // Save name of HDF5 file to be read to struct data
+  cross_sections_path_ = path_cross_sections;
+
   // Check if MGXS Library exists
-  if (!file_exists(settings::path_cross_sections)) {
+  if (!file_exists(cross_sections_path_)) {
     // Could not find MGXS Library file
-    fatal_error("Cross sections HDF5 file '" + settings::path_cross_sections +
+    fatal_error("Cross sections HDF5 file '" + cross_sections_path_ +
       "' does not exist.");
   }
   write_message("Reading cross sections HDF5 file...", 5);
 
   // Open file for reading
-  hid_t file_id = file_open(settings::path_cross_sections, 'r', true);
+  hid_t file_id = file_open(cross_sections_path_, 'r', true);
 
   ensure_exists(file_id, "energy_groups", true);
-  read_attribute(file_id, "energy_groups", data::num_energy_groups);
+  read_attribute(file_id, "energy_groups", num_energy_groups_);
 
   if (attribute_exists(file_id, "delayed_groups")) {
-    read_attribute(file_id, "delayed_groups", data::num_delayed_groups);
+    read_attribute(file_id, "delayed_groups", num_delayed_groups_);
   } else {
-    data::num_delayed_groups = 0;
+    num_delayed_groups_ = 0;
   }
 
   ensure_exists(file_id, "group structure", true);
-  read_attribute(file_id, "group structure", data::rev_energy_bins);
+  read_attribute(file_id, "group structure", rev_energy_bins_);
 
   // Reverse energy bins
-  std::copy(data::rev_energy_bins.crbegin(), data::rev_energy_bins.crend(),
-    std::back_inserter(data::energy_bins));
+  std::copy(rev_energy_bins_.crbegin(), rev_energy_bins_.crend(),
+    std::back_inserter(energy_bins_));
 
   // Create average energies
-  for (int i = 0; i < data::energy_bins.size() - 1; ++i) {
-    data::energy_bin_avg.push_back(0.5*(data::energy_bins[i] + data::energy_bins[i+1]));
+  for (int i = 0; i < energy_bins_.size() - 1; ++i) {
+    energy_bin_avg_.push_back(0.5*
+    (energy_bins_[i] + energy_bins_[i+1]));
   }
 
   // Add entries into libraries for MG data
-  auto names = group_names(file_id);
-  if (names.empty()) {
+  xs_names_ = group_names(file_id);
+  if (xs_names_.empty()) {
     fatal_error("At least one MGXS data set must be present in mgxs "
       "library file!");
   }
 
-  for (auto& name : names) {
+  // Close MGXS HDF5 file
+  file_close(file_id);
+}
+
+void put_mgxs_header_data_to_globals()
+{
+  // Get the minimum and maximum energies
+  int neutron = static_cast<int>(Particle::Type::neutron);
+  data::energy_min[neutron] = data::mg.energy_bins_.back();
+  data::energy_max[neutron] = data::mg.energy_bins_.front();
+
+  // Save available XS names to library list, so that when
+  // materials are read, the specified mgxs can be confirmed
+  // as present
+  for (auto& name : data::mg.xs_names_) {
     Library lib {};
     lib.type_ = Library::Type::neutron;
     lib.materials_.push_back(name);
     data::libraries.push_back(lib);
   }
+}
 
-  // Get the minimum and maximum energies
-  int neutron = static_cast<int>(Particle::Type::neutron);
-  data::energy_min[neutron] = data::energy_bins.back();
-  data::energy_max[neutron] = data::energy_bins.front();
+void set_mg_interface_nuclides_and_temps()
+{
+  // Get temperatures from global data
+  std::vector<std::vector<double>> nuc_temps(data::nuclide_map.size());
+  std::vector<std::vector<double>> dummy;
+  get_temperatures(nuc_temps, dummy);
 
-  // Close MGXS HDF5 file
-  file_close(file_id);
+  // Build vector of nuclide names which are to be read
+  std::vector<std::string> nuclide_names(data::nuclide_map.size());
+  for (const auto& kv : data::nuclide_map) {
+    nuclide_names[kv.second] = kv.first;
+  }
+
+  std::unordered_set<std::string> already_read;
+
+  // Loop over materials to find xs and temperature to be read
+  for (const auto& mat : model::materials) {
+    for (int i_nuc : mat->nuclide_) {
+      std::string& name = nuclide_names[i_nuc];
+
+      if (already_read.find(name) == already_read.end()) {
+        data::mg.xs_to_read_.push_back(name);
+        data::mg.xs_temps_to_read_.push_back(nuc_temps[i_nuc]);
+        already_read.insert(name);
+      }
+    }
+  }
+}
+
+void mark_fissionable_mgxs_materials()
+{
+  // Loop over all files
+  for (const auto& mat : model::materials) {
+    for (int i_nuc : mat->nuclide_) {
+      if (data::mg.nuclides_[i_nuc].fissionable) {
+        mat->fissionable_ = true;
+      }
+    }
+  }
 }
 
 //==============================================================================
@@ -251,7 +295,7 @@ void
 calculate_xs_c(int i_mat, int gin, double sqrtkT, Direction u,
      double& total_xs, double& abs_xs, double& nu_fiss_xs)
 {
-  data::macro_xs[i_mat].calculate_xs(gin - 1, sqrtkT, u, total_xs, abs_xs,
+  data::mg.macro_xs_[i_mat].calculate_xs(gin - 1, sqrtkT, u, total_xs, abs_xs,
        nu_fiss_xs);
 }
 
@@ -269,7 +313,7 @@ get_nuclide_xs(int index, int xstype, int gin, const int* gout,
   } else {
     gout_c_p = gout;
   }
-  return data::nuclides_MG[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg);
+  return data::mg.nuclides_[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg);
 }
 
 //==============================================================================
@@ -286,7 +330,7 @@ get_macro_xs(int index, int xstype, int gin, const int* gout,
   } else {
     gout_c_p = gout;
   }
-  return data::macro_xs[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg);
+  return data::mg.macro_xs_[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg);
 }
 
 //==============================================================================
@@ -301,7 +345,7 @@ get_name_c(int index, int name_len, char* name)
   std::strcpy(name, str.c_str());
 
   // Now get the data and copy to the C-string
-  str = data::nuclides_MG[index - 1].name;
+  str = data::mg.nuclides_[index - 1].name;
   std::strcpy(name, str.c_str());
 
   // Finally, remove the null terminator
@@ -313,7 +357,7 @@ get_name_c(int index, int name_len, char* name)
 double
 get_awr_c(int index)
 {
-  return data::nuclides_MG[index - 1].awr;
+  return data::mg.nuclides_[index - 1].awr;
 }
 
 } // namespace openmc
diff --git a/src/output.cpp b/src/output.cpp
index d3262cc4d5..899bce86e2 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -710,7 +710,7 @@ write_tallies()
               << data::nuclides[i_nuclide]->name_ << "\n";
           } else {
             tallies_out << std::string(indent+1, ' ')
-              << data::nuclides_MG[i_nuclide].name << "\n";
+              << data::mg.nuclides_[i_nuclide].name << "\n";
           }
         }
 
diff --git a/src/particle.cpp b/src/particle.cpp
index bf107a741c..030c8dc636 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -123,7 +123,7 @@ Particle::from_source(const Bank* src)
   } else {
     g_ = static_cast<int>(src->E);
     g_last_ = static_cast<int>(src->E);
-    E_ = data::energy_bin_avg[g_ - 1];
+    E_ = data::mg.energy_bin_avg_[g_ - 1];
   }
   E_last_ = E_;
 }
diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp
index b2017bbebe..8672ebf840 100644
--- a/src/particle_restart.cpp
+++ b/src/particle_restart.cpp
@@ -52,7 +52,7 @@ void read_particle_restart(Particle& p, int& previous_run_mode)
   // Set energy group and average energy in multi-group mode
   if (!settings::run_CE) {
     p.g_ = p.E_;
-    p.E_ = data::energy_bin_avg[p.g_ - 1];
+    p.E_ = data::mg.energy_bin_avg_[p.g_ - 1];
   }
 
   // Set particle last attributes
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index f99d15aac7..90e76c67a5 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -82,7 +82,7 @@ scatter(Particle* p)
   int gin = p->g_last_ - 1;
   int gout = p->g_ - 1;
   int i_mat = p->material_;
-  data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu_, p->wgt_);
+  data::mg.macro_xs_[i_mat].sample_scatter(gin, gout, p->mu_, p->wgt_);
 
   // Adjust return value for fortran indexing
   // TODO: Remove when no longer needed
@@ -92,7 +92,7 @@ scatter(Particle* p)
   p->u() = rotate_angle(p->u(), p->mu_, nullptr);
 
   // Update energy value for downstream compatability (in tallying)
-  p->E_ = data::energy_bin_avg[gout];
+  p->E_ = data::mg.energy_bin_avg_[gout];
 
   // Set event component
   p->event_ = EVENT_SCATTER;
@@ -148,7 +148,7 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
     // the energy in the fission bank
     int dg;
     int gout;
-    data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout);
+    data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_ - 1, dg, gout);
     site.E = gout + 1;
     site.delayed_group = dg + 1;
 
diff --git a/src/source.cpp b/src/source.cpp
index 0baa21b170..6b5cb8abac 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -311,9 +311,9 @@ Particle::Bank sample_external_source()
 
   // If running in MG, convert site % E to group
   if (!settings::run_CE) {
-    site.E = lower_bound_index(data::rev_energy_bins.begin(),
-      data::rev_energy_bins.end(), site.E);
-    site.E = data::num_energy_groups - site.E;
+    site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(),
+      data::mg.rev_energy_bins_.end(), site.E);
+    site.E = data::mg.num_energy_groups_ - site.E;
   }
 
   // Set the random number generator back to the tracking stream.
diff --git a/src/state_point.cpp b/src/state_point.cpp
index cd2c475d86..4c159c1678 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -208,7 +208,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
             if (settings::run_CE) {
               nuclides.push_back(data::nuclides[i_nuclide]->name_);
             } else {
-              nuclides.push_back(data::nuclides_MG[i_nuclide].name);
+              nuclides.push_back(data::mg.nuclides_[i_nuclide].name);
             }
           }
         }
diff --git a/src/summary.cpp b/src/summary.cpp
index d2721f97ef..1f5f7a097b 100644
--- a/src/summary.cpp
+++ b/src/summary.cpp
@@ -57,7 +57,7 @@ void write_nuclides(hid_t file)
       nuc_names.push_back(nuc->name_);
       awrs.push_back(nuc->awr_);
     } else {
-      const auto& nuc {data::nuclides_MG[i]};
+      const auto& nuc {data::mg.nuclides_[i]};
       if (nuc.awr != MACROSCOPIC_AWR) {
         nuc_names.push_back(nuc.name);
         awrs.push_back(nuc.awr);
diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp
index dde9b692ac..7737e4fac1 100644
--- a/src/tallies/filter_energy.cpp
+++ b/src/tallies/filter_energy.cpp
@@ -43,10 +43,10 @@ EnergyFilter::set_bins(gsl::span<const double> bins)
   // (after flipping for the different ordering of the library and tallying
   // systems).
   if (!settings::run_CE) {
-    if (n_bins_ == data::num_energy_groups) {
+    if (n_bins_ == data::mg.num_energy_groups_) {
       matches_transport_groups_ = true;
       for (gsl::index i = 0; i < n_bins_ + 1; ++i) {
-        if (data::rev_energy_bins[i] != bins_[i]) {
+        if (data::mg.rev_energy_bins_[i] != bins_[i]) {
           matches_transport_groups_ = false;
           break;
         }
@@ -61,9 +61,9 @@ const
 {
   if (p->g_ != F90_NONE && matches_transport_groups_) {
     if (estimator == ESTIMATOR_TRACKLENGTH) {
-      match.bins_.push_back(data::num_energy_groups - p->g_);
+      match.bins_.push_back(data::mg.num_energy_groups_ - p->g_);
     } else {
-      match.bins_.push_back(data::num_energy_groups - p->g_last_);
+      match.bins_.push_back(data::mg.num_energy_groups_ - p->g_last_);
     }
     match.weights_.push_back(1.0);
 
@@ -104,7 +104,7 @@ EnergyoutFilter::get_all_bins(const Particle* p, int estimator,
                               FilterMatch& match) const
 {
   if (p->g_ != F90_NONE && matches_transport_groups_) {
-    match.bins_.push_back(data::num_energy_groups - p->g_);
+    match.bins_.push_back(data::mg.num_energy_groups_ - p->g_);
     match.weights_.push_back(1.0);
 
   } else {
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 37310070fe..582e7dced9 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -361,7 +361,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
       if (settings::run_CE) {
         E_out = bank.E;
       } else {
-        E_out = data::energy_bin_avg[static_cast<int>(bank.E)];
+        E_out = data::mg.energy_bin_avg_[static_cast<int>(bank.E)];
       }
 
       // Set EnergyoutFilter bin index
@@ -1376,13 +1376,13 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
   // To significantly reduce de-referencing, point matxs to the macroscopic
   // Mgxs for the material of interest
-  data::macro_xs[p->material_].set_angle_index(p_u);
+  data::mg.macro_xs_[p->material_].set_angle_index(p_u);
 
   // Do same for nucxs, point it to the microscopic nuclide data of interest
   if (i_nuclide >= 0) {
     // And since we haven't calculated this temperature index yet, do so now
-    data::nuclides_MG[i_nuclide].set_temperature_index(p->sqrtkT_);
-    data::nuclides_MG[i_nuclide].set_angle_index(p_u);
+    data::mg.nuclides_[i_nuclide].set_temperature_index(p->sqrtkT_);
+    data::mg.nuclides_[i_nuclide].set_angle_index(p_u);
   }
 
   for (auto i = 0; i < tally.scores_.size(); ++i) {
@@ -1869,7 +1869,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
               // delayed-nu-fission xs to the absorption xs for all delayed
               // groups
               score = 0.;
-              for (auto d = 0; d < data::num_delayed_groups; ++d) {
+              for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) {
                 if (i_nuclide >= 0) {
                   score += p->wgt_absorb_ * flux
                     * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE,
@@ -1960,7 +1960,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           continue;
         } else {
           score = 0.;
-          for (auto d = 0; d < data::num_delayed_groups; ++d) {
+          for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) {
             if (i_nuclide >= 0) {
               score += atom_density * flux
                 * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE,
diff --git a/src/xsdata.cpp b/src/xsdata.cpp
index cee36a938b..d6b92958a1 100644
--- a/src/xsdata.cpp
+++ b/src/xsdata.cpp
@@ -24,11 +24,12 @@ namespace openmc {
 // XsData class methods
 //==============================================================================
 
-XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi)
+XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
+               size_t n_groups, size_t n_d_groups) :
+  n_g_(n_groups),
+  n_dg_(n_d_groups)
 {
   size_t n_ang = n_pol * n_azi;
-  size_t n_dg = data::num_delayed_groups;
-  size_t n_g = data::num_energy_groups;
 
   // check to make sure scatter format is OK before we allocate
   if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR &&
@@ -36,7 +37,7 @@ XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi)
     fatal_error("Invalid scatter_format!");
   }
   // allocate all [temperature][angle][in group] quantities
-  std::vector<size_t> shape {n_ang, n_g};
+  std::vector<size_t> shape {n_ang, n_g_};
   total = xt::zeros<double>(shape);
   absorption = xt::zeros<double>(shape);
   inverse_velocity = xt::zeros<double>(shape);
@@ -48,20 +49,20 @@ XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi)
   }
 
   // allocate decay_rate; [temperature][angle][delayed group]
-  shape[1] = n_dg;
+  shape[1] = n_dg_;
   decay_rate = xt::zeros<double>(shape);
 
   if (fissionable) {
-    shape = {n_ang, n_dg, n_g};
+    shape = {n_ang, n_dg_, n_g_};
     // 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]
-    shape = {n_ang, n_g, n_g};
+    shape = {n_ang, n_g_, n_g_};
     chi_prompt = xt::zeros<double>(shape);
 
     // chi_delayed; [temperature][angle][delay group][in group][out group]
-    shape = {n_ang, n_dg, n_g, n_g};
+    shape = {n_ang, n_dg_, n_g_, n_g_};
     chi_delayed = xt::zeros<double>(shape);
   }
 
@@ -127,11 +128,8 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
 {
   // Data is provided as nu-fission and chi with a beta for delayed info
 
-  size_t n_g = data::num_energy_groups;
-  size_t n_dg = data::num_delayed_groups;
-
   // Get chi
-  xt::xtensor<double, 2> temp_chi({n_ang, n_g}, 0.);
+  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
@@ -144,7 +142,7 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
                          xt::all());
 
   // Get nu-fission
-  xt::xtensor<double, 2> temp_nufiss({n_ang, n_g}, 0.);
+  xt::xtensor<double, 2> temp_nufiss({n_ang, n_g_}, 0.);
   read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true);
 
   // Get beta (strategy will depend upon the number of dimensions in beta)
@@ -154,7 +152,7 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   int ndim_target = 1;
   if (!is_isotropic) ndim_target += 2;
   if (beta_ndims == ndim_target) {
-    xt::xtensor<double, 2> temp_beta({n_ang, n_dg}, 0.);
+    xt::xtensor<double, 2> temp_beta({n_ang, n_dg_}, 0.);
     read_nd_vector(xsdata_grp, "beta", temp_beta, true);
 
     // Set prompt_nu_fission = (1. - beta_total)*nu_fission
@@ -165,7 +163,7 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
          xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) *
          xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all());
   } else if (beta_ndims == ndim_target + 1) {
-    xt::xtensor<double, 3> temp_beta({n_ang, n_dg, n_g}, 0.);
+    xt::xtensor<double, 3> temp_beta({n_ang, n_dg_, n_g_}, 0.);
     read_nd_vector(xsdata_grp, "beta", temp_beta, true);
 
     // Set prompt_nu_fission = (1. - beta_total)*nu_fission
@@ -182,18 +180,15 @@ 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
 
-  size_t n_g = data::num_energy_groups;
-  size_t n_dg = data::num_delayed_groups;
-
   // Get chi-prompt
-  xt::xtensor<double, 2> temp_chi_p({n_ang, n_g}, 0.);
+  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
   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.);
+  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
@@ -218,10 +213,8 @@ XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
   // No beta is provided and there is no prompt/delay distinction.
   // Therefore, the code only considers the data as prompt.
 
-  size_t n_g = data::num_energy_groups;
-
   // Get chi
-  xt::xtensor<double, 2> temp_chi({n_ang, n_g}, 0.);
+  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
@@ -241,11 +234,8 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_is
 {
   // Data is provided as nu-fission and chi with a beta for delayed info
 
-  size_t n_g = data::num_energy_groups;
-  size_t n_dg = data::num_delayed_groups;
-
   // Get nu-fission matrix
-  xt::xtensor<double, 3> temp_matrix({n_ang, n_g, n_g}, 0.);
+  xt::xtensor<double, 3> temp_matrix({n_ang, n_g_, n_g_}, 0.);
   read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true);
 
   // Get beta (strategy will depend upon the number of dimensions in beta)
@@ -255,7 +245,7 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_is
   int ndim_target = 1;
   if (!is_isotropic) ndim_target += 2;
   if (beta_ndims == ndim_target) {
-    xt::xtensor<double, 2> temp_beta({n_ang, n_dg}, 0.);
+    xt::xtensor<double, 2> temp_beta({n_ang, n_dg_}, 0.);
     read_nd_vector(xsdata_grp, "beta", temp_beta, true);
 
     xt::xtensor<double, 1> temp_beta_sum({n_ang}, 0.);
@@ -281,10 +271,10 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_is
          xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all());
 
   } else if (beta_ndims == ndim_target + 1) {
-    xt::xtensor<double, 3> temp_beta({n_ang, n_dg, n_g}, 0.);
+    xt::xtensor<double, 3> temp_beta({n_ang, n_dg_, n_g_}, 0.);
     read_nd_vector(xsdata_grp, "beta", temp_beta, true);
 
-    xt::xtensor<double, 2> temp_beta_sum({n_ang, n_g}, 0.);
+    xt::xtensor<double, 2> temp_beta_sum({n_ang, n_g_}, 0.);
     temp_beta_sum = xt::sum(temp_beta, {1});
 
     // prompt_nu_fission is the sum of this matrix over outgoing groups and
@@ -319,11 +309,8 @@ XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
 {
   // Data is provided separately as prompt + delayed nu-fission and chi
 
-  size_t n_g = data::num_energy_groups;
-  size_t n_dg = data::num_delayed_groups;
-
   // Get the prompt nu-fission matrix
-  xt::xtensor<double, 3> temp_matrix_p({n_ang, n_g, n_g}, 0.);
+  xt::xtensor<double, 3> temp_matrix_p({n_ang, n_g_, n_g_}, 0.);
   read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true);
 
   // prompt_nu_fission is the sum over outgoing groups
@@ -335,7 +322,7 @@ XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
        xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());
 
   // Get the delayed nu-fission matrix
-  xt::xtensor<double, 4> temp_matrix_d({n_ang, n_dg, n_g, n_g}, 0.);
+  xt::xtensor<double, 4> temp_matrix_d({n_ang, n_dg_, n_g_, n_g_}, 0.);
   read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_matrix_d, true);
 
   // delayed_nu_fission is the sum over outgoing groups
@@ -353,10 +340,8 @@ XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
   // No beta is provided and there is no prompt/delay distinction.
   // Therefore, the code only considers the data as prompt.
 
-  size_t n_g = data::num_energy_groups;
-
   // Get nu-fission matrix
-  xt::xtensor<double, 3> temp_matrix({n_ang, n_g, n_g}, 0.);
+  xt::xtensor<double, 3> temp_matrix({n_ang, n_g_, n_g_}, 0.);
   read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true);
 
   // prompt_nu_fission is the sum over outgoing groups
@@ -381,7 +366,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
   // as a nu-fission matrix or a set of chi and nu-fission vectors
   if (object_exists(xsdata_grp, "chi") ||
       object_exists(xsdata_grp, "chi-prompt")) {
-    if (data::num_delayed_groups == 0) {
+    if (n_dg_ == 0) {
       fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang);
     } else {
       if (object_exists(xsdata_grp, "beta")) {
@@ -391,7 +376,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
       }
     }
   } else {
-    if (data::num_delayed_groups == 0) {
+    if (n_dg_ == 0) {
       fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang);
     } else {
       if (object_exists(xsdata_grp, "beta")) {
@@ -403,7 +388,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
   }
 
   // Combine prompt_nu_fission and delayed_nu_fission into nu_fission
-  if (data::num_delayed_groups == 0) {
+  if (n_dg_ == 0) {
     nu_fission = prompt_nu_fission;
   } else {
     nu_fission = prompt_nu_fission + xt::sum(delayed_nu_fission, {1});
@@ -422,10 +407,9 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   hid_t scatt_grp = open_group(xsdata_grp, "scatter_data");
 
   // Get the outgoing group boundary indices
-  size_t n_g = data::num_energy_groups;
-  xt::xtensor<int, 2> gmin({n_ang, n_g}, 0.);
+  xt::xtensor<int, 2> gmin({n_ang, n_g_}, 0.);
   read_nd_vector(scatt_grp, "g_min", gmin, true);
-  xt::xtensor<int, 2> gmax({n_ang, n_g}, 0.);
+  xt::xtensor<int, 2> gmax({n_ang, n_g_}, 0.);
   read_nd_vector(scatt_grp, "g_max", gmax, true);
 
   // Make gmin and gmax start from 0 vice 1 as they do in the library
@@ -436,7 +420,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   // data.
   size_t length = order_data * xt::sum(gmax - gmin + 1)();
 
-  double_4dvec input_scatt(n_ang, double_3dvec(n_g));
+  double_4dvec input_scatt(n_ang, double_3dvec(n_g_));
   xt::xtensor<double, 1> temp_arr({length}, 0.);
   read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true);
 
@@ -453,7 +437,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   // scatt data
   size_t temp_idx = 0;
   for (size_t a = 0; a < n_ang; a++) {
-    for (size_t gin = 0; gin < n_g; gin++) {
+    for (size_t gin = 0; gin < n_g_; gin++) {
       input_scatt[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
       for (size_t i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) {
         input_scatt[a][gin][i_gout].resize(order_dim);
@@ -467,7 +451,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   }
 
   // Get multiplication matrix
-  double_3dvec temp_mult(n_ang, double_2dvec(n_g));
+  double_3dvec temp_mult(n_ang, double_2dvec(n_g_));
   if (object_exists(scatt_grp, "multiplicity_matrix")) {
     temp_arr.resize({length / order_data});
     read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr);
@@ -475,7 +459,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
     // convert the flat temp_arr to a jagged array for passing to scatt data
     size_t temp_idx = 0;
     for (size_t a = 0; a < n_ang; a++) {
-      for (size_t gin = 0; gin < n_g; gin++) {
+      for (size_t gin = 0; gin < n_g_; gin++) {
         temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
         for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
           temp_mult[a][gin][i_gout] = temp_arr[temp_idx++];
@@ -485,7 +469,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   } else {
     // Use a default: multiplicities are 1.0.
     for (size_t a = 0; a < n_ang; a++) {
-      for (size_t gin = 0; gin < n_g; gin++) {
+      for (size_t gin = 0; gin < n_g_; gin++) {
         temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
         for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
           temp_mult[a][gin][i_gout] = 1.;