diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml
index a3dae818e..5845898c8 100644
--- a/examples/xml/pincell/settings.xml
+++ b/examples/xml/pincell/settings.xml
@@ -3,9 +3,9 @@
 
   <!-- Define how many particles to run and for how many batches -->
   <run_mode>eigenvalue</run_mode>
-  <batches>20</batches>
+  <batches>100</batches>
   <inactive>10</inactive>
-  <particles>1000000</particles>
+  <particles>1000</particles>
 
   <!-- The starting source is a uniform distribution over the entire pin
        cell. Note that since this is effectively a 2D model, the z coordinates
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 00d37162b..7b0cd938d 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -330,6 +330,11 @@ public:
   std::vector<std::vector<Position>> tracks_;
 
   NuBank nu_bank_[15];
+
+  double tally_absorption_ {0.0};
+  double tally_collision_ {0.0};
+  double tally_tracklength_ {0.0};
+  double tally_leakage_ {0.0};
 };
 
 } // namespace openmc
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 71c83b5e8..5706b8ccb 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -167,8 +167,6 @@ extern double global_tally_absorption;
 extern double global_tally_collision;
 extern double global_tally_tracklength;
 extern double global_tally_leakage;
-#pragma omp threadprivate(global_tally_absorption, global_tally_collision, \
-  global_tally_tracklength, global_tally_leakage)
 
 //==============================================================================
 // Non-member functions
diff --git a/src/particle.cpp b/src/particle.cpp
index c4660b60a..e3f6db583 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -254,7 +254,7 @@ Particle::transport()
     // Score track-length estimate of k-eff
     if (settings::run_mode == RUN_MODE_EIGENVALUE &&
         type_ == Particle::Type::neutron) {
-      global_tally_tracklength += wgt_ * distance * macro_xs_.nu_fission;
+        tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission;
     }
 
     // Score flux derivative accumulators for differential tallies.
@@ -298,7 +298,7 @@ Particle::transport()
       // Score collision estimate of keff
       if (settings::run_mode == RUN_MODE_EIGENVALUE &&
           type_ == Particle::Type::neutron) {
-        global_tally_collision += wgt_ * macro_xs_.nu_fission
+        tally_collision_ += wgt_ * macro_xs_.nu_fission
           / macro_xs_.total;
       }
 
@@ -398,6 +398,16 @@ Particle::transport()
     write_particle_track(*this);
     finalize_particle_track(*this);
   }
+
+  // Contribute tally reduction variables to global accumulator
+  #pragma omp atomic
+  global_tally_absorption += tally_absorption_;
+  #pragma omp atomic
+  global_tally_collision += tally_collision_;
+  #pragma omp atomic
+  global_tally_tracklength += tally_tracklength_;
+  #pragma omp atomic
+  global_tally_leakage += tally_leakage_;
 }
 
 void
@@ -430,7 +440,7 @@ Particle::cross_surface()
     }
 
     // Score to global leakage tally
-    global_tally_leakage += wgt_;
+    tally_leakage_ += wgt_;
 
     // Display message
     if (settings::verbosity >= 10 || simulation::trace) {
diff --git a/src/physics.cpp b/src/physics.cpp
index ac73f9c85..c639c5907 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -573,7 +573,7 @@ void absorption(Particle* p, int i_nuclide)
 
     // Score implicit absorption estimate of keff
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      global_tally_absorption += p->wgt_absorb_ * p->neutron_xs_[
+      p->tally_absorption_ += p->wgt_absorb_ * p->neutron_xs_[
         i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
     }
   } else {
@@ -582,7 +582,7 @@ void absorption(Particle* p, int i_nuclide)
         prn(p->current_seed()) * p->neutron_xs_[i_nuclide].total) {
       // Score absorption estimate of keff
       if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-        global_tally_absorption += p->wgt_ * p->neutron_xs_[
+        p->tally_absorption_ += p->wgt_ * p->neutron_xs_[
           i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
       }
 
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index bc25abc6f..f2bd8cc42 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -177,13 +177,11 @@ absorption(Particle* p)
     p->wgt_last_ = p->wgt_;
 
     // Score implicit absorpion estimate of keff
-    #pragma omp atomic
-    global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission /
+    p->tally_absorption_ += p->wgt_absorb_ * p->macro_xs_.nu_fission /
         p->macro_xs_.absorption;
   } else {
     if (p->macro_xs_.absorption > prn(p->current_seed()) * p->macro_xs_.total) {
-      #pragma omp atomic
-      global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission /
+      p->tally_absorption_ += p->wgt_ * p->macro_xs_.nu_fission /
            p->macro_xs_.absorption;
       p->alive_ = false;
       p->event_ = EVENT_ABSORB;
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 1f86d0850..f0dd4ad07 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -354,7 +354,7 @@ void process_advance_particle_events()
     // 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;
+      p->tally_tracklength_ += p->wgt_ * distance * p->macro_xs_.nu_fission;
     }
 
     // Score flux derivative accumulators for differential tallies.
@@ -422,7 +422,7 @@ void process_collision_events()
     // 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->tally_collision_ += p->wgt_ * p->macro_xs_.nu_fission
         / p->macro_xs_.total;
     }
 
@@ -713,13 +713,17 @@ void transport()
 			}
 		}
 
-    // Finish particle track output.
+    // Finish particle track output and contribute to global tally variables
     for (int i = 0; i < n_particles; i++) {
       Particle& p = particles[i];
       if (p.write_track_) {
         write_particle_track(p);
         finalize_particle_track(p);
       }
+      global_tally_absorption +=  p.tally_absorption_;
+      global_tally_collision +=   p.tally_collision_;
+      global_tally_tracklength += p.tally_tracklength_;
+      global_tally_leakage +=     p.tally_leakage_;
     }
 
 		remaining_work -= n_particles;
@@ -904,7 +908,7 @@ int openmc_next_batch(int* status)
 
     // ====================================================================
     // LOOP OVER PARTICLES
-
+    /*
     #pragma omp parallel for schedule(runtime)
     for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
       // grab source particle from bank
@@ -914,8 +918,9 @@ int openmc_next_batch(int* status)
       // transport particle
       p.transport();
     }
+    */
 
-    //transport();
+    transport();
 
     // Accumulate time for transport
     simulation::time_transport.stop();
@@ -1218,6 +1223,12 @@ void initialize_history(Particle* p, int64_t index_source)
   p->from_source(&simulation::source_bank[index_source - 1]);
   p->current_work_ = index_source;
 
+  // Initialize global tally variables
+  p->tally_absorption_ =  0.0;
+  p->tally_collision_ =   0.0;
+  p->tally_tracklength_ = 0.0;
+  p->tally_leakage_ =     0.0;
+
   // set identifier for particle
   p->id_ = simulation::work_index[mpi::rank] + index_source;