diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml
index 5845898c8..7b0a985f1 100644
--- a/examples/xml/pincell/settings.xml
+++ b/examples/xml/pincell/settings.xml
@@ -3,8 +3,8 @@
 
   <!-- Define how many particles to run and for how many batches -->
   <run_mode>eigenvalue</run_mode>
-  <batches>100</batches>
-  <inactive>10</inactive>
+  <batches>10</batches>
+  <inactive>5</inactive>
   <particles>1000</particles>
 
   <!-- The starting source is a uniform distribution over the entire pin
diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h
index e386348c4..a6e1e4063 100644
--- a/include/openmc/geometry.h
+++ b/include/openmc/geometry.h
@@ -28,12 +28,14 @@ extern std::vector<int64_t> overlap_check_count;
 // Information about nearest boundary crossing
 //==============================================================================
 
+/*
 struct BoundaryInfo {
   double distance {INFINITY};   //!< distance to nearest boundary
   int surface_index {0}; //!< if boundary is surface, index in surfaces vector
   int coord_level;   //!< coordinate level after crossing boundary
   std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
 };
+*/
 
 //==============================================================================
 //! Check two distances by coincidence tolerance
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 6ac8aabae..760e2ad8a 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -129,8 +129,21 @@ struct MacroXS {
   double incoherent;      //!< macroscopic incoherent xs
   double photoelectric;   //!< macroscopic photoelectric xs
   double pair_production; //!< macroscopic pair production xs
+  int i_grid;
 };
 
+//==============================================================================
+// Information about nearest boundary crossing
+//==============================================================================
+
+struct BoundaryInfo {
+  double distance {INFINITY};   //!< distance to nearest boundary
+  int surface_index {0}; //!< if boundary is surface, index in surfaces vector
+  int coord_level;   //!< coordinate level after crossing boundary
+  std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
+};
+
+
 //============================================================================
 //! State of a particle being transported through geometry
 //============================================================================
@@ -273,6 +286,7 @@ public:
 
   // Indices for various arrays
   int surface_ {0};             //!< index for surface particle is on
+  BoundaryInfo boundary_;
   int cell_born_ {-1};      //!< index for cell particle was born in
   int material_ {-1};       //!< index for current material
   int material_last_ {-1};  //!< index for last material
@@ -294,4 +308,14 @@ public:
 
 } // namespace openmc
 
+extern template class std::vector<openmc::Particle::Bank>;
+
+namespace openmc{
+
+extern std::vector<Particle> particle_bank;
+#pragma omp threadprivate(particle_bank)
+
+}
+
+
 #endif // OPENMC_PARTICLE_H
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index 8efac4c7f..7981286be 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -16,6 +16,10 @@ constexpr int STATUS_EXIT_NORMAL {0};
 constexpr int STATUS_EXIT_MAX_BATCH {1};
 constexpr int STATUS_EXIT_ON_TRIGGER {2};
 
+extern Particle::Bank * shared_fission_bank;
+extern int shared_fission_bank_length;
+extern int shared_fission_bank_max;
+
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -38,7 +42,9 @@ extern "C" int restart_batch;   //!< batch at which a restart job resumed
 extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
 extern "C" int total_gen;        //!< total number of generations simulated
 extern double total_weight;  //!< Total source weight in a batch
-extern int64_t work_per_rank;         //!< number of particles per MPI rank
+extern int64_t thread_work_index;
+extern int64_t work_per_rank;       //!< number of particles per MPI rank
+extern int64_t work_per_thread;    //!< number of particles on a given tread
 
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;
@@ -49,7 +55,8 @@ extern std::vector<int64_t> work_index;
 // Threadprivate variables
 extern "C" bool trace;     //!< flag to show debug information
 
-#pragma omp threadprivate(current_work, trace)
+
+#pragma omp threadprivate(current_work, work_per_thread, trace)
 
 } // namespace simulation
 
diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index 59263afc8..f946d0602 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -85,6 +85,7 @@ double AngleDistribution::sample(double E, uint64_t * prn_seeds, int stream) con
   // Sample between the ith and (i+1)th bin
   if (r > prn(prn_seeds, stream)) ++i;
 
+  //assert(i >= 0 );
   // Sample i-th distribution
   double mu = distribution_[i]->sample(prn_seeds, stream);
 
diff --git a/src/material.cpp b/src/material.cpp
index a5bc3dc13..16b9275a0 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -751,6 +751,7 @@ void Material::calculate_xs(Particle& p) const
 
 void Material::calculate_neutron_xs(Particle& p) const
 {
+	/*
   // Find energy index on energy grid
   int neutron = static_cast<int>(Particle::Type::neutron);
   int i_grid = std::log(p.E_/data::energy_min[neutron])/simulation::log_spacing;
@@ -760,9 +761,11 @@ void Material::calculate_neutron_xs(Particle& p) const
 
   // Initialize position in i_sab_nuclides
   int j = 0;
+  */
 
   // Add contribution from each nuclide in material
   for (int i = 0; i < nuclide_.size(); ++i) {
+	  /*
     // ======================================================================
     // CHECK FOR S(A,B) TABLE
 
@@ -792,18 +795,21 @@ void Material::calculate_neutron_xs(Particle& p) const
 
     // ======================================================================
     // CALCULATE MICROSCOPIC CROSS SECTION
+	*/
 
     // Determine microscopic cross sections for this nuclide
     int i_nuclide = nuclide_[i];
 
     // Calculate microscopic cross section for this nuclide
     const auto& micro {p.neutron_xs_[i_nuclide]};
+	/*
     if (p.E_ != micro.last_E
         || p.sqrtkT_ != micro.last_sqrtkT
         || i_sab != micro.index_sab
         || sab_frac != micro.sab_frac) {
       data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p);
     }
+	*/
 
     // ======================================================================
     // ADD TO MACROSCOPIC CROSS SECTION
diff --git a/src/output.cpp b/src/output.cpp
index 899bce86e..29d60830e 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -67,7 +67,15 @@ void title()
     "                     ###############     %%%%%%%%%%%%%%%%\n" <<
     "                       ############     %%%%%%%%%%%%%%%\n" <<
     "                          ########     %%%%%%%%%%%%%%\n" <<
-    "                                      %%%%%%%%%%%\n\n";
+    "                                      %%%%%%%%%%%\n" <<
+    "                              )                             (        (\n" << 
+    "                           ( /(   *   )        (     (      )\\ )     )\\ )\n" <<  
+    "         (    (   (   (    )\\())` )  /(      ( )\\    )\\    (()/( (  (()/(\n" << 
+    "         )\\   )\\  )\\  )\\  ((_)\\  ( )(_))___  )((_)((((_)(   /(_)))\\  /(_))\n" << 
+    "        ((_) ((_)((_)((_)  _((_)(_(_())|___|((_)_  )\\ _ )\\ (_)) ((_)(_))_\n" <<  
+    "        | __|\\ \\ / / | __|| \\| ||_   _|      | _ ) (_)_\\(_)/ __|| __||   \\\n" << 
+    "        | _|  \\ V /  | _| | .` |  | |        | _ \\  / _ \\  \\__ \\| _| | |) |\n" << 
+    "        |___|  \\_/   |___||_|\\_|  |_|        |___/ /_/ \\_\\ |___/|___||___/\n\n";
 
   // Write version information
   std::cout <<
diff --git a/src/particle.cpp b/src/particle.cpp
index 98ee8920b..12334f149 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -27,8 +27,14 @@
 #include "openmc/tallies/tally_scoring.h"
 #include "openmc/track_output.h"
 
+// Explicit template instantiation definition
+template class std::vector<openmc::Particle>;
+
+
 namespace openmc {
 
+std::vector<Particle> particle_bank;
+
 //==============================================================================
 // LocalCoord implementation
 //==============================================================================
@@ -132,14 +138,15 @@ void
 Particle::transport()
 {
   // Display message if high verbosity or trace is on
-  if (settings::verbosity >= 9 || simulation::trace) {
-     write_message("Simulating Particle " + std::to_string(id_));
-  }
+  //if (settings::verbosity >= 9 || simulation::trace) {
+    // write_message("Simulating Particle " + std::to_string(id_));
+  //}
 
   // Initialize number of events to zero
   int n_event = 0;
 
   // Add paricle's starting weight to count for normalizing tallies later
+  /*
   #pragma omp atomic
   simulation::total_weight += wgt_;
 
@@ -153,6 +160,7 @@ Particle::transport()
 
   // Every particle starts with no accumulated flux derivative.
   if (!model::active_tallies.empty()) zero_flux_derivs();
+  */
 
   while (true) {
     // Set the random number stream
diff --git a/src/physics.cpp b/src/physics.cpp
index 28fc4341e..84927c27e 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -79,6 +79,7 @@ void collision(Particle* p)
 
 void sample_neutron_reaction(Particle* p)
 {
+	//std::cout << "particle energy e = " << p->E_ << std::endl;
   // Sample a nuclide within the material
   int i_nuclide = sample_nuclide(p);
 
@@ -127,9 +128,12 @@ void sample_neutron_reaction(Particle* p)
   }
   if (!p->alive_) return;
 
+	//std::cout << "before scatter particle energy e = " << p->E_ << " with nuclide = " << i_nuclide <<  std::endl;
   // Sample a scattering reaction and determine the secondary energy of the
   // exiting neutron
   scatter(p, i_nuclide);
+	
+  //std::cout << "after scatter particle energy e = " << p->E_ << std::endl;
 
   // Advance URR seed stream 'N' times after energy changes
   if (p->E_ != p->E_last_) {
@@ -171,6 +175,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 
   p->fission_ = true;
   for (int i = 0; i < nu; ++i) {
+	  /*
     // Create new bank site and get reference to last element
     bank.emplace_back();
     auto& site {bank.back()};
@@ -179,12 +184,22 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
     site.r = p->r();
     site.particle = Particle::Type::neutron;
     site.wgt = 1. / weight;
+	*/
+	  int idx;
+	  #pragma omp atomic capture
+	  idx = shared_fission_bank_length++;
+	  Particle::Bank * site = shared_fission_bank + idx;
+	  site->r = p->r();
+	  site->particle = Particle::Type::neutron;
+	  site->wgt = 1. / weight;
 
     // Sample delayed group and angle/energy for fission reaction
-    sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->prn_seeds_, p->stream_);
+
+    sample_fission_neutron(i_nuclide, rx, p->E_, site, p->prn_seeds_, p->stream_);
 
     // Set the delayed group on the particle as well
-    p->delayed_group_ = site.delayed_group;
+    //p->delayed_group_ = site.delayed_group;
+    p->delayed_group_ = site->delayed_group;
 
     // Increment the number of neutrons born delayed
     if (p->delayed_group_ > 0) {
@@ -598,6 +613,7 @@ void scatter(Particle* p, int i_nuclide)
     // Determine temperature
     double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp];
 
+	//std::cout << "we are doing an elastic scatter" << std::endl;
     // Perform collision physics for elastic scattering
     elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p);
 
@@ -610,6 +626,7 @@ void scatter(Particle* p, int i_nuclide)
     // =======================================================================
     // S(A,B) SCATTERING
 
+	//std::cout << "we are doing an SAB scatter" << std::endl;
     sab_scatter(i_nuclide, micro.index_sab, p);
 
     p->event_mt_ = ELASTIC;
@@ -617,6 +634,7 @@ void scatter(Particle* p, int i_nuclide)
   }
 
   if (!sampled) {
+	//std::cout << "we are doing an Inelastic scatter" << std::endl;
     // =======================================================================
     // INELASTIC SCATTERING
 
@@ -643,7 +661,9 @@ void scatter(Particle* p, int i_nuclide)
 
     // Perform collision physics for inelastic scattering
     const auto& rx {nuc->reactions_[i]};
+	//std::cout << "Energy before Inelastic scatter = " << p->E_ << std::endl;
     inelastic_scatter(nuc.get(), rx.get(), p);
+	//std::cout << "Energy after Inelastic scatter = " << p->E_ << std::endl;
     p->event_mt_ = rx->mt_;
   }
 
@@ -1056,22 +1076,29 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
   // angle and outgoing energy from CM to LAB
   if (rx->scatter_in_cm_) {
     double E_cm = E;
+	//std::cout << "E_cm = " << E_cm << std::endl;
 
     // determine outgoing energy in lab
     double A = nuc->awr_;
+//	std::cout << "A = " << nuc->awr_ << std::endl;
     E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm))
           / ((A + 1.0)*(A + 1.0));
+	// here's where it goes wrong
+	//std::cout << "E = " << E << std::endl;
 
     // determine outgoing angle in lab
     mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E);
+//	std::cout << "mu = " << mu << std::endl;
   }
 
   // Because of floating-point roundoff, it may be possible for mu to be
   // outside of the range [-1,1). In these cases, we just set mu to exactly -1
   // or 1
   if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
+//	std::cout << "mu = " << mu << std::endl;
 
   // Set outgoing energy and scattering angle
+//	std::cout << "FINAL E = " << E << " mu = " << mu << std::endl;
   p->E_ = E;
   p->mu_ = mu;
 
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 86f20a93d..b2af14a3f 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -2,23 +2,32 @@
 
 #include "openmc/bank.h"
 #include "openmc/capi.h"
+#include "openmc/cell.h"
 #include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
+#include "openmc/geometry.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
+#include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
 #include "openmc/particle.h"
 #include "openmc/photon.h"
+#include "openmc/physics.h"
+ #include "openmc/physics_mg.h"
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
 #include "openmc/source.h"
 #include "openmc/state_point.h"
+#include "openmc/thermal.h"
 #include "openmc/timer.h"
+#include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/tallies/tally_scoring.h"
 #include "openmc/tallies/trigger.h"
+#include "openmc/track_output.h"
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -28,6 +37,604 @@
 #include <algorithm>
 #include <string>
 
+namespace openmc {
+
+	/*
+extern std::vector<Particle*> calculate_fuel_xs_queue;
+extern std::vector<Particle*> calculate_nonfuel_xs_queue;
+extern std::vector<Particle*> advance_particle_queue;
+extern std::vector<Particle*> surface_crossing_queue;
+extern std::vector<Particle*> collision_queue;
+#pragma omp threadprivate(calculate_fuel_xs_queue, calculate_nonfuel_xs_queue, advance_particle_queue, surface_crossing_queue, collision_queue)
+
+std::vector<Particle*> calculate_fuel_xs_queue;
+std::vector<Particle*> calculate_nonfuel_xs_queue;
+std::vector<Particle*> advance_particle_queue;
+std::vector<Particle*> surface_crossing_queue;
+std::vector<Particle*> collision_queue;
+*/
+
+int * calculate_fuel_xs_queue;
+int * calculate_nonfuel_xs_queue;
+int * advance_particle_queue;
+int * surface_crossing_queue;
+int * collision_queue;
+Particle * particles;
+
+int calculate_fuel_xs_queue_length    = 0;
+int calculate_nonfuel_xs_queue_length = 0;
+int advance_particle_queue_length     = 0;
+int surface_crossing_queue_length     = 0;
+int collision_queue_length            = 0;
+
+const int MAX_PARTICLES_IN_FLIGHT = 100;
+
+void init_event_queues(void)
+{
+	int n_particles = MAX_PARTICLES_IN_FLIGHT;
+	calculate_fuel_xs_queue =    new int[n_particles];
+	calculate_nonfuel_xs_queue = new int[n_particles];
+	advance_particle_queue =     new int[n_particles];
+	surface_crossing_queue =     new int[n_particles];
+	collision_queue =            new int[n_particles];
+	particles =                  new Particle[n_particles];
+}
+
+void free_event_queues(void)
+{
+	delete calculate_fuel_xs_queue;
+	delete calculate_nonfuel_xs_queue;
+	delete advance_particle_queue;
+	delete surface_crossing_queue;
+	delete collision_queue;
+	delete[] particles;
+}
+
+constexpr size_t MAX_PARTICLES_PER_THREAD {100};
+
+Particle::Bank * shared_fission_bank;
+int shared_fission_bank_length = 0;
+int shared_fission_bank_max;
+
+void init_shared_fission_bank(int max)
+{
+	shared_fission_bank_max = max;
+	shared_fission_bank = new Particle::Bank[max];
+}
+
+void free_shared_fission_bank(void)
+{
+	delete[] shared_fission_bank;
+	shared_fission_bank_length = 0;
+}
+
+// TODO: What is going on here?
+void revive_particle_from_secondary(Particle* p)
+{
+  p->from_source(&simulation::secondary_bank.back());
+  simulation::secondary_bank.pop_back();
+  // n_event = 0;
+
+  // Enter new particle in particle track file
+  if (p->write_track_) add_particle_track();
+}
+
+void dispatch_xs_event(int i)
+{
+	Particle * p = particles + i;
+	  int idx;
+  if (p->material_ == MATERIAL_VOID) {
+
+	#pragma omp atomic capture
+	idx = calculate_nonfuel_xs_queue_length++;
+	//std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl;
+    calculate_nonfuel_xs_queue[idx] = i;
+  } else {
+    if (model::materials[p->material_]->fissionable_) {
+	#pragma omp atomic capture
+	idx = calculate_fuel_xs_queue_length++;
+	//std::cout << "Dispatching particle to Fuel XS queue idx = " << idx << std::endl;
+      calculate_fuel_xs_queue[idx] = i;
+    } else {
+	#pragma omp atomic capture
+	idx = calculate_nonfuel_xs_queue_length++;
+	//std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl;
+      calculate_nonfuel_xs_queue[idx] = i;
+    }
+  }
+}
+
+void process_calculate_xs_events(int * queue, int n)
+{
+  // Save last_ members, find grid index
+  for (int i = 0; i < n; i++) {
+	  Particle *p = particles + queue[i]; 
+	  //std::cout << "particle offset = " << queue[i] << std::endl;
+    // Set the random number stream
+	// TODO: Move RNG seeds to particle storage
+    if (p->type_ == Particle::Type::neutron) {
+      prn_set_stream(STREAM_TRACKING);
+    } else {
+      prn_set_stream(STREAM_PHOTON);
+    }
+
+    // Store pre-collision particle properties
+    p->wgt_last_ = p->wgt_;
+    p->E_last_ = p->E_;
+    p->u_last_ = p->u();
+    p->r_last_ = p->r();
+
+    // If the cell hasn't been determined based on the particle's location,
+    // initiate a search for the current cell. This generally happens at the
+    // beginning of the history and again for any secondary particles
+    if (p->coord_[p->n_coord_ - 1].cell == C_NONE) {
+      if (!find_cell(p, false)) {
+        p->mark_as_lost("Could not find the cell containing particle "
+          + std::to_string(p->id_));
+        return;
+      }
+
+      // set birth cell attribute
+      if (p->cell_born_ == C_NONE) p->cell_born_ = p->coord_[p->n_coord_ - 1].cell;
+    }
+
+    // Write particle track.
+    if (p->write_track_) write_particle_track(*p);
+
+    if (settings::check_overlaps) check_cell_overlap(p);
+
+    if (settings::run_CE) {
+      if (p->material_ == p->material_last_ && p->sqrtkT_ != p->sqrtkT_last_) {
+        // Remove particle from queue
+      }
+    }
+
+    // Find energy index on energy grid
+    // TODO: Calculate this separately?
+    int neutron = static_cast<int>(Particle::Type::neutron);
+    p->macro_xs_.i_grid = std::log(p->E_/data::energy_min[neutron]) / simulation::log_spacing;
+  }
+
+  // Calculate nuclide micros
+  for (int i = 0; i < data::nuclides.size(); ++i) {
+	  // loop over particles
+    for (int j = 0; j < n; j++) {
+		//Particle * p = particles + queue[i];
+		Particle * p = particles + queue[j];
+      if (p->material_ == MATERIAL_VOID) continue;
+
+      // If material doesn't have this nuclide, skip it
+      const auto& mat {model::materials[p->material_]};
+      if (mat->mat_nuclide_index_[i] == -1) continue;
+
+      // ======================================================================
+      // CHECK FOR S(A,B) TABLE
+
+      // Check if this nuclide matches one of the S(a,b) tables specified.
+      // This relies on thermal_tables_ being sorted by .index_nuclide
+      int i_sab = C_NONE;
+      double sab_frac = 0.0;
+      for (const auto& sab : mat->thermal_tables_) {
+        if (i == sab.index_nuclide) {
+          // Get index in sab_tables
+          i_sab = sab.index_table;
+          sab_frac = sab.fraction;
+
+          // If particle energy is greater than the highest energy for the
+          // S(a,b) table, then don't use the S(a,b) table
+          //if (p->E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE;
+          if (p->E_ > data::thermal_scatt[i_sab]->energy_max_) i_sab = C_NONE;
+        }
+      }
+
+      // ======================================================================
+      // CALCULATE MICROSCOPIC CROSS SECTION
+
+      // Calculate microscopic cross section for this nuclide
+      const auto& micro {p->neutron_xs_[i]};
+      if (p->E_ != micro.last_E
+          || p->sqrtkT_ != micro.last_sqrtkT
+          || i_sab != micro.index_sab
+          || sab_frac != micro.sab_frac) {
+        data::nuclides[i]->calculate_xs(i_sab, p->macro_xs_.i_grid, sab_frac, *p);
+      }
+    }
+  }
+
+  for (int i = 0; i < n; i++) {
+	  Particle * p = particles + queue[i];
+    // Calculate microscopic and macroscopic cross sections
+    if (p->material_ != MATERIAL_VOID) {
+	  // Only works for CE, no MG support
+      //if (settings::run_CE) {
+        // If the material is the same as the last material and the
+        // temperature hasn't changed, we don't need to lookup cross
+        // sections again.
+        model::materials[p->material_]->calculate_xs(*p);
+      }
+		//else {
+      //}
+     else {
+      p->macro_xs_.total      = 0.0;
+      p->macro_xs_.absorption = 0.0;
+      p->macro_xs_.fission    = 0.0;
+      p->macro_xs_.nu_fission = 0.0;
+    }
+
+	 int idx;
+	 #pragma omp atomic capture
+	 idx = advance_particle_queue_length++;
+    advance_particle_queue[idx] = queue[i];
+  }
+}
+
+void process_advance_particle_events()
+{
+  //for (auto& p : advance_particle_queue) {
+  for (int i = 0; i < advance_particle_queue_length; i++) {
+	  Particle * p = particles + advance_particle_queue[i];
+    simulation::trace == (p->id_ == 0);
+
+    // Sample a distance to collision
+    double d_collision;
+    if (p->type_ == Particle::Type::electron ||
+        p->type_ == Particle::Type::positron) {
+      d_collision = 0.0;
+    } else if (p->macro_xs_.total == 0.0) {
+      d_collision = INFINITY;
+    } else {
+      d_collision = -std::log(prn()) / p->macro_xs_.total;
+    }
+
+    // -------------- break here? -------------------
+
+    // Find the distance to the nearest boundary
+    p->boundary_ = distance_to_boundary(p);
+
+    // Select smaller of the two distances
+    double distance;
+	int idx;
+    if (p->boundary_.distance < d_collision) {
+		#pragma omp atomic capture
+		idx = surface_crossing_queue_length++;
+      surface_crossing_queue[idx] = advance_particle_queue[i];
+      distance = p->boundary_.distance;
+    } else {
+		#pragma omp atomic capture
+		idx = collision_queue_length++;
+		collision_queue[idx] = advance_particle_queue[i];
+      distance = d_collision;
+    }
+
+    // -------------- break here? -------------------
+
+    // Advance particle
+    for (int j = 0; j < p->n_coord_; ++j) {
+      p->coord_[j].r += distance * p->coord_[j].u;
+    }
+
+    // -------------- break here? -------------------
+
+    // Score track-length tallies
+    if (!model::active_tracklength_tallies.empty()) {
+      score_tracklength_tally(p, distance);
+    }
+
+    // Score track-length estimate of k-eff
+    if (settings::run_mode == RUN_MODE_EIGENVALUE &&
+        p->type_ == Particle::Type::neutron) {
+      global_tally_tracklength += p->wgt_ * distance * p->macro_xs_.nu_fission;
+    }
+
+    // Score flux derivative accumulators for differential tallies.
+    if (!model::active_tallies.empty()) {
+      score_track_derivative(p, distance);
+    }
+  }
+
+  advance_particle_queue_length = 0;
+}
+
+void process_surface_crossing_events()
+{
+  //for (auto& p : surface_crossing_queue) {
+  for (int i = 0; i < surface_crossing_queue_length; i++) {
+	  Particle * p = particles + surface_crossing_queue[i];
+    // Set surface that particle is on and adjust coordinate levels
+    p->surface_ = p->boundary_.surface_index;
+    p->n_coord_ = p->boundary_.coord_level;
+
+    // Saving previous cell data
+    for (int j = 0; j < p->n_coord_; ++j) {
+      p->cell_last_[j] = p->coord_[j].cell;
+    }
+    p->n_coord_last_ = p->n_coord_;
+
+    if (p->boundary_.lattice_translation[0] != 0 ||
+        p->boundary_.lattice_translation[1] != 0 ||
+        p->boundary_.lattice_translation[2] != 0) {
+      // Particle crosses lattice boundary
+      cross_lattice(p, p->boundary_);
+      p->event_ = EVENT_LATTICE;
+    } else {
+      // Particle crosses surface
+      p->cross_surface();
+      p->event_ = EVENT_SURFACE;
+    }
+    // Score cell to cell partial currents
+    if (!model::active_surface_tallies.empty()) {
+      score_surface_tally(p, model::active_surface_tallies);
+    }
+
+	// TODO: Add back in support for secondary particles
+	/*
+    if (!p->alive_ && !simulation::secondary_bank.empty()) {
+      revive_particle_from_secondary(p);
+    }
+	*/
+
+    if (p->alive_)
+	{
+		dispatch_xs_event(surface_crossing_queue[i]);
+	}
+  }
+
+  surface_crossing_queue_length = 0;
+}
+
+void process_collision_events()
+{
+  //for (auto& p : collision_queue) {
+  for (int i = 0; i < collision_queue_length; i++) {
+	  Particle * p = particles + collision_queue[i];
+	  //std::cout << "Beginning collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
+    // Score collision estimate of keff
+    if (settings::run_mode == RUN_MODE_EIGENVALUE &&
+        p->type_ == Particle::Type::neutron) {
+      global_tally_collision += p->wgt_ * p->macro_xs_.nu_fission
+        / p->macro_xs_.total;
+    }
+
+    // Score surface current tallies -- this has to be done before the collision
+    // since the direction of the particle will change and we need to use the
+    // pre-collision direction to figure out what mesh surfaces were crossed
+
+    if (!model::active_meshsurf_tallies.empty())
+      score_surface_tally(p, model::active_meshsurf_tallies);
+	    
+	//std::cout << "After surface tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
+
+    // Clear surface component
+    p->surface_ = 0;
+
+    if (settings::run_CE) {
+      collision(p);
+    } else {
+      collision_mg(p);
+    }
+
+	//std::cout << "After collision() of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
+
+    // Score collision estimator tallies -- this is done after a collision
+    // has occurred rather than before because we need information on the
+    // outgoing energy for any tallies with an outgoing energy filter
+    if (!model::active_collision_tallies.empty()) score_collision_tally(p);
+    if (!model::active_analog_tallies.empty()) {
+      if (settings::run_CE) {
+        score_analog_tally_ce(p);
+      } else {
+        score_analog_tally_mg(p);
+      }
+    }
+	//std::cout << "After analog tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
+
+    // Reset banked weight during collision
+    p->n_bank_ = 0;
+    p->wgt_bank_ = 0.0;
+    for (int& v : p->n_delayed_bank_) v = 0;
+
+    // Reset fission logical
+    p->fission_ = false;
+
+    // Save coordinates for tallying purposes
+    p->r_last_current_ = p->r();
+
+    // Set last material to none since cross sections will need to be
+    // re-evaluated
+    p->material_last_ = C_NONE;
+
+    // Set all directions to base level -- right now, after a collision, only
+    // the base level directions are changed
+    for (int j = 0; j < p->n_coord_ - 1; ++j) {
+      if (p->coord_[j + 1].rotated) {
+        // If next level is rotated, apply rotation matrix
+        const auto& m {model::cells[p->coord_[j].cell]->rotation_};
+        const auto& u {p->coord_[j].u};
+        p->coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z;
+        p->coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z;
+        p->coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z;
+      } else {
+        // Otherwise, copy this level's direction
+        p->coord_[j+1].u = p->coord_[j].u;
+      }
+    }
+
+    // Score flux derivative accumulators for differential tallies.
+    if (!model::active_tallies.empty()) score_collision_derivative(p);
+
+	// TODO: Add back in secondary particle support
+	/*
+    if (!p->alive_ && !simulation::secondary_bank.empty()) {
+      revive_particle_from_secondary(p);
+    }
+	*/
+
+    if (p->alive_)
+	{
+		dispatch_xs_event(collision_queue[i]);
+	    //std::cout << "Ended collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
+		assert(std::isfinite(p->E_) );
+	}
+  }
+
+  collision_queue_length = 0;
+}
+
+void check_energies(void)
+{
+	int * Q;
+	int n;
+
+
+	Q = calculate_fuel_xs_queue;
+	n = calculate_fuel_xs_queue_length;
+	for( int i = 0; i < n; i++ )
+	{
+		if( !std::isfinite(particles[Q[i]].E_ ) )
+		{
+			std::cout << "NAN energy particle found at index xs FUEL " << Q[i] << std::endl;
+			assert(0);
+		}
+	}
+	Q = calculate_nonfuel_xs_queue;
+	n = calculate_nonfuel_xs_queue_length ;
+	for( int i = 0; i < n; i++ )
+	{
+		if( !std::isfinite(particles[Q[i]].E_ ) )
+		{
+			std::cout << "NAN energy particle found at index xs Non fuel " << Q[i] << std::endl;
+			assert(0);
+		}
+	}
+	Q = advance_particle_queue;
+	n = advance_particle_queue_length     ;
+	for( int i = 0; i < n; i++ )
+	{
+		if( !std::isfinite(particles[Q[i]].E_ ) )
+		{
+			std::cout << "NAN energy particle found at index advance particle " << Q[i] << std::endl;
+			assert(0);
+		}
+	}
+	Q = surface_crossing_queue;
+	n = surface_crossing_queue_length     ;
+	for( int i = 0; i < n; i++ )
+	{
+		if( !std::isfinite(particles[Q[i]].E_ ) )
+		{
+			std::cout << "NAN energy particle found at index surface crossing " << Q[i] << std::endl;
+			assert(0);
+		}
+	}
+	Q = collision_queue;
+	n = collision_queue_length            ;
+	for( int i = 0; i < n; i++ )
+	{
+		if( !std::isfinite(particles[Q[i]].E_ ) )
+		{
+			std::cout << "NAN energy particle found at index collision " << Q[i] << std::endl;
+			assert(0);
+		}
+	}
+}
+
+void transport()
+{
+	int remaining_work = simulation::work_per_rank;
+	int source_offset = 0;
+
+
+
+	// Subiterations to complete sets of particles
+	while (remaining_work > 0) {
+
+		// TODO: Do this only once per PI
+		init_event_queues();
+
+		// Figure out work for this subiteration
+		int n_particles = MAX_PARTICLES_IN_FLIGHT;
+		if( n_particles > remaining_work)
+			n_particles = remaining_work;
+
+		//std::cout << "Initializing particle histories..." << std::endl;
+		// Initialize all histories
+		// TODO: Parallelize
+		for (int i = 0; i < n_particles; i++) {
+			initialize_history(particles + i, source_offset + i);
+		}
+
+		//std::cout << "Enqueing particles for XS Lookups..." << std::endl;
+		// Add all particles to advance particle queue
+		// TODO: Parallelize
+		for (int i = 0; i < n_particles; i++) {
+			dispatch_xs_event(i);
+		}
+
+		int event_kernel_executions = 0;
+		while (true) {
+			event_kernel_executions++;
+			/*
+			std::cout << "Fuel XS Lookups = " << calculate_fuel_xs_queue_length << std::endl;
+			std::cout << "Non Fuel XS Lookups = " << calculate_nonfuel_xs_queue_length << std::endl;
+			std::cout << "Advance Particles = " << advance_particle_queue_length << std::endl;
+			std::cout << "Surface Crossings = " << surface_crossing_queue_length << std::endl;
+			std::cout << "Collisions = " << collision_queue_length << std::endl;
+			*/
+			int max = std::max({calculate_fuel_xs_queue_length, calculate_nonfuel_xs_queue_length, advance_particle_queue_length, surface_crossing_queue_length, collision_queue_length});
+			check_energies();
+			if (max == 0) {
+				break;
+			} else if (max == calculate_fuel_xs_queue_length) {
+				//std::cout << "pre fuel XS check..." << std::endl;
+				//check_energies(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
+				//std::cout << "Performing Fuel XS Lookups..." << std::endl;
+				process_calculate_xs_events(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
+				calculate_fuel_xs_queue_length = 0;
+			} else if (max == calculate_nonfuel_xs_queue_length) {
+				//std::cout << "pre non fuel XS check..." << std::endl;
+				//check_energies(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
+			//	std::cout << "Performing Non Fuel XS Lookups..." << std::endl;
+				process_calculate_xs_events(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
+				calculate_nonfuel_xs_queue_length = 0;
+			} else if (max == advance_particle_queue_length) {
+				//std::cout << "pre advancing check..." << std::endl;
+				//check_energies(advance_particle_queue, advance_particle_queue_length);
+				//std::cout << "Advancing Particles..." << std::endl;
+				process_advance_particle_events();
+			} else if (max == surface_crossing_queue_length) {
+				//std::cout << "pre surface crossing check..." << std::endl;
+				//check_energies(surface_crossing_queue, surface_crossing_queue_length);
+				//std::cout << "Surface Crossings..." << std::endl;
+				process_surface_crossing_events();
+			} else if (max == collision_queue_length) {
+				//std::cout << "pre Colliding check..." << std::endl;
+				//check_energies(collision_queue, collision_queue_length);
+				//std::cout << "Colliding..." << std::endl;
+				process_collision_events();
+			}
+		}
+
+		remaining_work -= n_particles;
+		source_offset += n_particles;
+
+		// Should all be zero
+		/*
+		calculate_fuel_xs_queue_length    = 0;
+		calculate_nonfuel_xs_queue_length = 0;
+		advance_particle_queue_length     = 0;
+		surface_crossing_queue_length     = 0;
+		collision_queue_length            = 0;
+		*/
+		
+		// TODO: Do this only once per PI
+		free_event_queues();
+		std::cout << "Event kernels retired: " << event_kernel_executions << std::endl;
+	}
+	//shared_fission_bank_length = 0;
+}
+
+} // namespace openmc
+
 //==============================================================================
 // C API functions
 //==============================================================================
@@ -42,10 +649,12 @@ int openmc_run()
 
   int err = 0;
   int status = 0;
+
   while (status == 0 && err == 0) {
     err = openmc_next_batch(&status);
   }
 
+
   openmc_simulation_finalize();
   return err;
 }
@@ -56,9 +665,16 @@ int openmc_simulation_init()
 
   // Skip if simulation has already been initialized
   if (simulation::initialized) return 0;
+	
 
   // Determine how much work each process should do
   calculate_work();
+  
+  // We assume that the number of fission sites will be at most
+	// 10x the number of neutrons run on this rank
+	int shared_fission_bank_size = 10 * simulation::work_per_rank;
+	// Initializes the shared fission bank
+	init_shared_fission_bank(shared_fission_bank_size);
 
   // Allocate array for matching filter bins
   #pragma omp parallel
@@ -71,9 +687,11 @@ int openmc_simulation_init()
   allocate_banks();
 
   // Allocate tally results arrays if they're not allocated yet
+	  std::cout << "Trying to Initialize Results..." << std::endl;
   for (auto& t : model::tallies) {
     t->init_results();
   }
+	  std::cout << "Success in Initializing Results!" << std::endl;
 
   // Set up material nuclide index mapping
   for (auto& mat : model::materials) {
@@ -118,6 +736,8 @@ int openmc_simulation_finalize()
 
   // Skip if simulation was never run
   if (!simulation::initialized) return 0;
+  
+  free_shared_fission_bank();
 
   // Stop active batch timer and start finalization timer
   simulation::time_active.stop();
@@ -188,6 +808,14 @@ int openmc_next_batch(int* status)
     // ====================================================================
     // LOOP OVER PARTICLES
 
+	simulation::current_work = 1;
+
+    //#pragma omp parallel
+    {
+      transport();
+    }
+
+/*
     #pragma omp parallel for schedule(runtime)
     for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
       simulation::current_work = i_work;
@@ -199,6 +827,7 @@ int openmc_next_batch(int* status)
       // transport particle
       p.transport();
     }
+*/
 
     // Accumulate time for transport
     simulation::time_transport.stop();
@@ -255,7 +884,9 @@ int restart_batch;
 bool satisfy_triggers {false};
 int total_gen {0};
 double total_weight;
+int64_t thread_work_index;
 int64_t work_per_rank;
+int64_t work_per_thread;
 
 const RegularMesh* entropy_mesh {nullptr};
 const RegularMesh* ufs_mesh {nullptr};
@@ -439,7 +1070,15 @@ void finalize_generation()
   if (settings::run_mode == RUN_MODE_EIGENVALUE) {
 #ifdef _OPENMP
     // Join the fission bank from each thread into one global fission bank
-    join_bank_from_threads();
+    //join_bank_from_threads();
+	
+	  // We need to move all the stuff from the shared_fission_bank into the real one.
+		//std::vector<Particle::Bank> shared_fission_bank_vector(shared_fission_bank, shared_fission_bank + shared_fission_bank_length);
+      //simulation::fission_bank = shared_fission_bank_vector;
+	  std::cout << "Fission bank length = " << shared_fission_bank_length << std::endl;
+	  for( int i = 0; i < shared_fission_bank_length; i++ )
+		  simulation::fission_bank.push_back(shared_fission_bank[i]);
+	  shared_fission_bank_length = 0;
 #endif
 
     // Distribute fission bank across processors evenly
@@ -467,8 +1106,10 @@ void finalize_generation()
 
 void initialize_history(Particle* p, int64_t index_source)
 {
+	//assert(index_source - 1 >= 0 );
   // set defaults
-  p->from_source(&simulation::source_bank[index_source - 1]);
+  //p->from_source(&simulation::source_bank[index_source - 1]);
+  p->from_source(&simulation::source_bank[index_source]);
 
   // set identifier for particle
   p->id_ = simulation::work_index[mpi::rank] + index_source;
@@ -498,6 +1139,29 @@ void initialize_history(Particle* p, int64_t index_source)
       }
     }
   }
+
+// Display message if high verbosity or trace is on
+   if (settings::verbosity >= 9 || simulation::trace) {
+      write_message("Simulating Particle " + std::to_string(p->id_));
+   }
+
+    // // Initialize number of events to zero
+   // int n_event = 0;
+
+    // Add paricle's starting weight to count for normalizing tallies later
+   #pragma omp atomic
+   simulation::total_weight += p->wgt_;
+
+    // Force calculation of cross-sections by setting last energy to zero
+   if (settings::run_CE) {
+     for (auto& micro : p->neutron_xs_) micro.last_E = 0.0;
+   }
+
+    // Prepare to write out particle track.
+   if (p->write_track_) add_particle_track();
+
+    // Every particle starts with no accumulated flux derivative.
+   if (!model::active_tallies.empty()) zero_flux_derivs();
 }
 
 int overall_generation()
@@ -528,6 +1192,32 @@ void calculate_work()
     i_bank += work_i;
     simulation::work_index[i + 1] = i_bank;
   }
+	#ifdef _OPENMP
+   // Determine work per thread
+   int remainder_thread = simulation::work_per_rank % omp_get_max_threads();
+
+    #pragma omp parallel
+   {
+     simulation::work_per_thread = simulation::work_per_rank / omp_get_num_threads();
+     if (omp_get_thread_num() < remainder_thread) {
+       ++simulation::work_per_thread;
+     }
+   }
+
+    int64_t work_i = 0;
+   #pragma omp parallel for ordered
+   for (int i = 0; i < omp_get_num_threads(); ++i) {
+     #pragma omp ordered
+     {
+       simulation::thread_work_index = work_i;
+       work_i += simulation::work_per_thread;
+	   /*
+       std::cout << "Thread " << omp_get_thread_num() << ": " << simulation::work_per_thread
+         << std::endl;
+		 */
+     }
+   }
+ #endif
 }
 
 #ifdef OPENMC_MPI
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 852918a0c..bfd91df1d 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -804,6 +804,7 @@ Tally::init_triggers(pugi::xml_node node)
 void Tally::init_results()
 {
   int n_scores = scores_.size() * nuclides_.size();
+  std::cout << "TRYING TO ALLOCATE n filter bins, n scores =  " << n_filter_bins_ << " " <<  n_scores << std::endl;
   results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
 }