diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h
index cce7ca164f..690cb5ae87 100644
--- a/include/openmc/hdf5_interface.h
+++ b/include/openmc/hdf5_interface.h
@@ -52,7 +52,7 @@ void write_string(hid_t group_id, const char* name, const std::string& buffer,
 
 std::vector<hsize_t> attribute_shape(hid_t obj_id, const char* name);
 std::vector<std::string> dataset_names(hid_t group_id);
-void ensure_exists(hid_t group_id, const char* name);
+void ensure_exists(hid_t obj_id, const char* name, bool attribute=false);
 std::vector<std::string> group_names(hid_t group_id);
 std::vector<hsize_t> object_shape(hid_t obj_id);
 std::string object_name(hid_t obj_id);
diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h
index e0addf5bc0..c813755111 100644
--- a/include/openmc/mgxs_interface.h
+++ b/include/openmc/mgxs_interface.h
@@ -7,6 +7,7 @@
 #include "hdf5_interface.h"
 #include "mgxs.h"
 
+#include <vector>
 
 namespace openmc {
 
@@ -17,9 +18,9 @@ namespace openmc {
 extern std::vector<Mgxs> nuclides_MG;
 extern std::vector<Mgxs> macro_xs;
 extern "C" int num_energy_groups;
-
-//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also
-// bring over energy_bin_avg, energy_bins, etc, as vectors)
+extern std::vector<double> energy_bins;
+extern std::vector<double> energy_bin_avg;
+extern std::vector<double> rev_energy_bins;
 
 //==============================================================================
 // Mgxs data loading interface methods
@@ -39,6 +40,8 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
      int n_temps, const double temps[], const double atom_densities[],
      double tolerance, int& method);
 
+extern "C" void read_mg_cross_sections_header_c(hid_t file_id);
+
 //==============================================================================
 // Mgxs tracking/transport/tallying interface methods
 //==============================================================================
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 298b2878a0..6f2d58832b 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -146,9 +146,7 @@ extern "C" {
     //! site may have been produced from an external source, from fission, or
     //! simply as a secondary particle.
     //! \param src Source site data
-    //! \param run_CE Whether continuous-energy data is being used
-    //! \param energy_bin_avg An array of energy group bin averages
-    void from_source(const Bank* src, bool run_CE, const double* energy_bin_avg);
+    void from_source(const Bank* src);
 
     //! mark a particle as lost and create a particle restart file
     //! \param message A warning message to display
@@ -174,8 +172,7 @@ extern "C" {
   void particle_create_secondary(Particle* p, const double* uvw, double E,
                                  int type, bool run_CE);
   void particle_initialize(Particle* p);
-  void particle_from_source(Particle* p, const Bank* src, bool run_CE,
-                            const double* energy_bin_avg);
+  void particle_from_source(Particle* p, const Bank* src);
   void particle_mark_as_lost(Particle* p, const char* message);
   void particle_write_restart(Particle* p);
 
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 6565d3e460..2a3111fe67 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -7,6 +7,7 @@
 #include <array>
 #include <cstdint>
 #include <string>
+#include <vector>
 
 #include "pugixml.hpp"
 
@@ -71,7 +72,6 @@ extern "C" int legendre_to_tabular_points; //!< number of points to convert Lege
 extern "C" int max_order;                //!< Maximum Legendre order for multigroup data
 extern "C" int n_log_bins;               //!< number of bins for logarithmic energy grid
 extern "C" int n_max_batches;            //!< Maximum number of batches
-
 extern "C" int res_scat_method;          //!< resonance upscattering method
 extern "C" double res_scat_energy_min;   //!< Min energy in [eV] for res. upscattering
 extern "C" double res_scat_energy_max;   //!< Max energy in [eV] for res. upscattering
@@ -83,6 +83,7 @@ extern "C" double temperature_range[2];  //!< Min/max T in [K] over which to loa
 extern "C" int trace_batch;              //!< Batch to trace particle on
 extern "C" int trace_gen;                //!< Generation to trace particle on
 extern "C" int64_t trace_particle;       //!< Particle ID to enable trace on
+extern std::vector<std::array<int, 3>> track_identifiers; //!< Particle numbers for writing tracks
 extern "C" int trigger_batch_interval;   //!< Batch interval for triggers
 extern "C" int verbosity;                //!< How verbose to make output
 extern "C" double weight_cutoff;         //!< Weight cutoff for Russian roulette
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index a6279f285c..71b2e7e741 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -4,6 +4,8 @@
 #ifndef OPENMC_SIMULATION_H
 #define OPENMC_SIMULATION_H
 
+#include "openmc/particle.h"
+
 #include <cstdint>
 #include <vector>
 
@@ -62,6 +64,8 @@ extern "C" void initialize_batch();
 //! Initialize a fission generation
 extern "C" void initialize_generation();
 
+extern "C" void initialize_history(Particle* p, int64_t index_source);
+
 //! Finalize a fission generation
 extern "C" void finalize_generation();
 
diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp
index d72f7d424d..4b0aedd0ae 100644
--- a/src/hdf5_interface.cpp
+++ b/src/hdf5_interface.cpp
@@ -148,13 +148,22 @@ dataset_typesize(hid_t dset)
 
 
 void
-ensure_exists(hid_t group_id, const char* name)
+ensure_exists(hid_t obj_id, const char* name, bool attribute)
 {
-  if (!object_exists(group_id, name)) {
-    std::stringstream err_msg;
-    err_msg << "Object \"" << name << "\" does not exist in group "
-            << object_name(group_id);
-    fatal_error(err_msg);
+  if (attribute) {
+    if (!attribute_exists(obj_id, name)) {
+      std::stringstream err_msg;
+      err_msg << "Attribute \"" << name << "\" does not exist in object "
+              << object_name(obj_id);
+      fatal_error(err_msg);
+    }
+  } else {
+    if (!object_exists(obj_id, name)) {
+      std::stringstream err_msg;
+      err_msg << "Object \"" << name << "\" does not exist in object "
+              << object_name(obj_id);
+      fatal_error(err_msg);
+    }
   }
 }
 
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index d86afd3cf1..2a8ded7339 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -211,7 +211,6 @@ contains
     integer :: i
     integer :: n
     integer, allocatable :: temp_int_array(:)
-    integer :: n_tracks
     type(XMLNode) :: root
     type(XMLNode) :: node_sp
     type(XMLNode) :: node_res_scat
@@ -221,25 +220,6 @@ contains
     ! Get proper XMLNode type given pointer
     root % ptr = root_ptr
 
-    ! Particle tracks
-    if (check_for_node(root, "track")) then
-      ! Make sure that there are three values per particle
-      n_tracks = node_word_count(root, "track")
-      if (mod(n_tracks, 3) /= 0) then
-        call fatal_error("Number of integers specified in 'track' is not &
-             &divisible by 3.  Please provide 3 integers per particle to be &
-             &tracked.")
-      end if
-
-      ! Allocate space and get list of tracks
-      allocate(temp_int_array(n_tracks))
-      call get_node_array(root, "track", temp_int_array)
-
-      ! Reshape into track_identifiers
-      allocate(track_identifiers(3, n_tracks/3))
-      track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
-    end if
-
     ! Check if the user has specified to write state points
     if (check_for_node(root, "state_point")) then
 
@@ -2719,6 +2699,13 @@ contains
     integer(HID_T) :: file_id
     character(len=MAX_WORD_LEN), allocatable :: names(:)
 
+    interface
+      subroutine read_mg_cross_sections_header_c(file_id) bind(C)
+        import HID_T
+        integer(HID_T), value :: file_id
+      end subroutine
+    end interface
+
     ! Check if MGXS Library exists
     inquire(FILE=path_cross_sections, EXIST=file_exists)
     if (.not. file_exists) then
@@ -2766,6 +2753,9 @@ contains
       energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1))
     end do
 
+    ! Set up energy bins on C++ side
+    call read_mg_cross_sections_header_c(file_id)
+
     ! Get the minimum and maximum energies
     energy_min(NEUTRON) = energy_bins(num_energy_groups + 1)
     energy_max(NEUTRON) = energy_bins(1)
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index 12c73fab4f..1b6e85ffb2 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -26,7 +26,6 @@ namespace openmc {
 std::vector<Mgxs> nuclides_MG;
 std::vector<Mgxs> macro_xs;
 
-
 //==============================================================================
 // Mgxs base-class methods
 //==============================================================================
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index c592d28a28..5cec8811f1 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -8,6 +8,14 @@
 
 namespace openmc {
 
+//==============================================================================
+// Global variable definitions
+//==============================================================================
+
+std::vector<double> energy_bins;
+std::vector<double> energy_bin_avg;
+std::vector<double> rev_energy_bins;
+
 //==============================================================================
 // Mgxs data loading interface methods
 //==============================================================================
@@ -83,6 +91,26 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
   }
 }
 
+//==============================================================================
+
+void read_mg_cross_sections_header_c(hid_t file_id)
+{
+  ensure_exists(file_id, "energy_groups", true);
+  read_attribute(file_id, "energy_groups", num_energy_groups);
+
+  ensure_exists(file_id, "group structure", true);
+  read_attribute(file_id, "group structure", energy_bins);
+
+  // Create reverse energy bins
+  std::copy(energy_bins.crbegin(), energy_bins.crend(),
+    std::back_inserter(rev_energy_bins));
+
+  // Create average energies
+  for (int i = 0; i < energy_bins.size() - 1; ++i) {
+    energy_bin_avg.push_back(0.5*(energy_bins[i] + energy_bins[i+1]));
+  }
+}
+
 //==============================================================================
 // Mgxs tracking/transport/tallying interface methods
 //==============================================================================
diff --git a/src/particle.cpp b/src/particle.cpp
index d6dfa969c1..c8a5a1b567 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -7,6 +7,7 @@
 #include "openmc/constants.h"
 #include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/mgxs_interface.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 
@@ -92,7 +93,7 @@ Particle::initialize()
 }
 
 void
-Particle::from_source(const Bank* src, bool run_CE, const double* energy_bin_avg)
+Particle::from_source(const Bank* src)
 {
   // set defaults
   initialize();
@@ -106,7 +107,7 @@ Particle::from_source(const Bank* src, bool run_CE, const double* energy_bin_avg
   std::copy(src->xyz, src->xyz + 3, last_xyz_current);
   std::copy(src->xyz, src->xyz + 3, last_xyz);
   std::copy(src->uvw, src->uvw + 3, last_uvw);
-  if (run_CE) {
+  if (settings::run_CE) {
     E = src->E;
     g = 0;
   } else {
@@ -212,10 +213,9 @@ void particle_create_secondary(Particle* p, const double* uvw, double E,
   p->create_secondary(uvw, E, type, run_CE);
 }
 void particle_initialize(Particle* p) { p->initialize(); }
-void particle_from_source(Particle* p, const Bank* src, bool run_CE,
-                          const double* energy_bin_avg)
+void particle_from_source(Particle* p, const Bank* src)
 {
-  p->from_source(src, run_CE, energy_bin_avg);
+  p->from_source(src);
 }
 void particle_mark_as_lost(Particle* p, const char* message)
 {
diff --git a/src/particle_header.F90 b/src/particle_header.F90
index df5deff8c8..db2317f4b2 100644
--- a/src/particle_header.F90
+++ b/src/particle_header.F90
@@ -129,12 +129,10 @@ module particle_header
       type(Particle), intent(inout) :: p
     end subroutine particle_initialize
 
-    subroutine particle_from_source(p, src, run_CE, energy_bin_avg) bind(C)
-      import Particle, Bank, C_BOOL, C_DOUBLE
+    subroutine particle_from_source(p, src) bind(C)
+      import Particle, Bank, C_DOUBLE
       type(Particle), intent(inout) :: p
       type(Bank),     intent(in)    :: src
-      logical(C_BOOL), value        :: run_CE
-      real(C_DOUBLE), intent(in)    :: energy_bin_avg(*)
     end subroutine particle_from_source
 
     subroutine particle_mark_as_lost_c(p, message) bind(C, name='particle_mark_as_lost')
diff --git a/src/settings.F90 b/src/settings.F90
index 5929bdd939..d0353342b4 100644
--- a/src/settings.F90
+++ b/src/settings.F90
@@ -99,7 +99,6 @@ module settings
 
   ! Particle tracks
   logical(C_BOOL), bind(C) :: write_all_tracks
-  integer, allocatable :: track_identifiers(:,:)
 
   ! Particle restart run
   logical(C_BOOL), bind(C) :: particle_restart_run
@@ -151,7 +150,6 @@ contains
 
   subroutine free_memory_settings()
     if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
-    if (allocated(track_identifiers)) deallocate(track_identifiers)
 
     call statepoint_batch % clear()
     call sourcepoint_batch % clear()
diff --git a/src/settings.cpp b/src/settings.cpp
index 556fbd24db..3442c9a5a2 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -94,6 +94,7 @@ double temperature_range[2] {0.0, 0.0};
 int trace_batch;
 int trace_gen;
 int64_t trace_particle;
+std::vector<std::array<int, 3>> track_identifiers;
 int trigger_batch_interval {1};
 int verbosity {7};
 double weight_cutoff {0.25};
@@ -495,7 +496,8 @@ void read_settings_xml()
   // Particle tracks
   if (check_for_node(root, "track")) {
     // Get values and make sure there are three per particle
-    auto temp = get_node_array<int64_t>(root, "track");
+    auto temp = get_node_array<int>(root, "track");
+    std::cout << "track size " << temp.size() << "\n";
     if (temp.size() % 3 != 0) {
       fatal_error("Number of integers specified in 'track' is not "
         "divisible by 3.  Please provide 3 integers per particle to be "
@@ -503,8 +505,11 @@ void read_settings_xml()
     }
 
     // Reshape into track_identifiers
-    //allocate(track_identifiers(3, n_tracks/3))
-    //track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
+    int n_tracks = temp.size() / 3;
+    for (int i = 0; i < n_tracks; ++i) {
+      settings::track_identifiers.push_back({temp[3*i], temp[3*i + 1],
+        temp[3*i + 2]});
+    }
   }
 
   // Read meshes
diff --git a/src/simulation.F90 b/src/simulation.F90
index 2bcdf68ab6..18fadc9a5f 100644
--- a/src/simulation.F90
+++ b/src/simulation.F90
@@ -55,6 +55,12 @@ module simulation
     subroutine finalize_generation() bind(C)
     end subroutine finalize_generation
 
+    subroutine initialize_history(p, index_source) bind(C)
+      import Particle, C_INT64_T
+      type(Particle), intent(inout) :: p
+      integer(C_INT64_T), value :: index_source
+    end subroutine
+
     function sample_external_source() result(site) bind(C)
       import Bank
       type(Bank) :: site
@@ -135,51 +141,6 @@ contains
 
   end function openmc_next_batch
 
-!===============================================================================
-! INITIALIZE_HISTORY
-!===============================================================================
-
-  subroutine initialize_history(p, index_source)
-
-    type(Particle), intent(inout) :: p
-    integer(8),     intent(in)    :: index_source
-
-    integer(8) :: particle_seed  ! unique index for particle
-    integer :: i
-
-    ! set defaults
-    call particle_from_source(p, source_bank(index_source), run_CE, &
-         energy_bin_avg)
-
-    ! set identifier for particle
-    p % id = work_index(rank) + index_source
-
-    ! set random number seed
-    particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id
-    call set_particle_seed(particle_seed)
-
-    ! set particle trace
-    trace = .false.
-    if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
-         p % id == trace_particle) trace = .true.
-
-    ! Set particle track.
-    p % write_track = .false.
-    if (write_all_tracks) then
-      p % write_track = .true.
-    else if (allocated(track_identifiers)) then
-      do i=1, size(track_identifiers(1,:))
-        if (current_batch == track_identifiers(1,i) .and. &
-             &current_gen == track_identifiers(2,i) .and. &
-             &p % id == track_identifiers(3,i)) then
-          p % write_track = .true.
-          exit
-        end if
-      end do
-    end if
-
-  end subroutine initialize_history
-
 !===============================================================================
 ! FINALIZE_BATCH handles synchronization and accumulation of tallies,
 ! calculation of Shannon entropy, getting single-batch estimate of keff, and
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 60c5737401..c4f1f03872 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -5,6 +5,7 @@
 #include "openmc/error.h"
 #include "openmc/message_passing.h"
 #include "openmc/output.h"
+#include "openmc/random_lcg.h"
 #include "openmc/settings.h"
 #include "openmc/source.h"
 #include "openmc/timer.h"
@@ -303,6 +304,46 @@ void finalize_generation()
   }
 }
 
+void initialize_history(Particle* p, int64_t index_source)
+{
+  // Get pointer to source bank
+  Bank* source_bank;
+  int64_t n;
+  openmc_source_bank(&source_bank, &n);
+
+  // set defaults
+  p->from_source(&source_bank[index_source - 1]);
+
+  // set identifier for particle
+  p->id = simulation::work_index[mpi::rank] + index_source;
+
+  // set random number seed
+  int64_t particle_seed = (simulation::total_gen + overall_generation() - 1)
+    * settings::n_particles + p->id;
+  set_particle_seed(particle_seed);
+
+  // set particle trace
+  simulation::trace = false;
+  if (simulation::current_batch == settings::trace_batch &&
+      simulation::current_gen == settings::trace_gen &&
+      p->id == settings::trace_particle) simulation::trace = true;
+
+  // Set particle track.
+  p->write_track = false;
+  if (settings::write_all_tracks) {
+    p->write_track = true;
+  } else if (settings::track_identifiers.size() > 0) {
+    for (const auto& t : settings::track_identifiers) {
+      if (simulation::current_batch == t[0] &&
+          simulation::current_gen == t[1] &&
+          p->id == t[2]) {
+        p->write_track = true;
+        break;
+      }
+    }
+  }
+}
+
 int overall_generation()
 {
   using namespace simulation;
diff --git a/src/tracking.F90 b/src/tracking.F90
index dc0c008d61..d14de4a461 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -287,8 +287,7 @@ contains
       ! Check for secondary particles if this particle is dead
       if (.not. p % alive) then
         if (p % n_secondary > 0) then
-          call particle_from_source(p, p % secondary_bank(p % n_secondary), &
-                                    run_CE, energy_bin_avg)
+          call particle_from_source(p, p % secondary_bank(p % n_secondary))
           p % n_secondary = p % n_secondary - 1
           n_event = 0