diff --git a/include/openmc/bank.h b/include/openmc/bank.h
index b7fab9ac66..4597369c78 100644
--- a/include/openmc/bank.h
+++ b/include/openmc/bank.h
@@ -20,15 +20,13 @@ namespace openmc {
 
 namespace simulation {
 
-extern "C" int64_t n_bank;
-
 extern std::vector<Particle::Bank> source_bank;
 extern std::vector<Particle::Bank> fission_bank;
 #ifdef _OPENMP
 extern std::vector<Particle::Bank> master_fission_bank;
 #endif
 
-#pragma omp threadprivate(fission_bank, n_bank)
+#pragma omp threadprivate(fission_bank)
 
 } // namespace simulation
 
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 3fc7db7816..f53bba224d 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -96,14 +96,11 @@ public:
 
   //! Count number of bank sites in each mesh bin / energy bin
   //
-  //! \param[in] n Number of bank sites
   //! \param[in] bank Array of bank sites
-  //! \param[in] n_energy Number of energies
-  //! \param[in] energies Array of energies
   //! \param[out] Whether any bank sites are outside the mesh
   //! \return Array indicating number of sites in each mesh/energy bin
-  xt::xarray<double> count_sites(int64_t n, const Particle::Bank* bank,
-    int n_energy, const double* energies, bool* outside) const;
+  xt::xarray<double> count_sites(const std::vector<Particle::Bank>& bank,
+    bool* outside) const;
 
   int id_ {-1};  //!< User-specified ID
   int n_dimension_; //!< Number of dimensions
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index ef1e076bb6..d52c67ea55 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -285,8 +285,7 @@ public:
   bool write_track_ {false};
 
   // Secondary particles created
-  int64_t n_secondary_ {};
-  Bank secondary_bank_[MAX_SECONDARY];
+  std::vector<Bank> secondary_bank_;
 };
 
 } // namespace openmc
diff --git a/include/openmc/physics.h b/include/openmc/physics.h
index 5c8d64f5f9..107380eb84 100644
--- a/include/openmc/physics.h
+++ b/include/openmc/physics.h
@@ -7,6 +7,8 @@
 #include "openmc/position.h"
 #include "openmc/reaction.h"
 
+#include <vector>
+
 namespace openmc {
 
 //==============================================================================
@@ -47,7 +49,7 @@ int sample_nuclide(const Particle* p);
 //! Determine the average total, prompt, and delayed neutrons produced from
 //! fission and creates appropriate bank sites.
 void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
-  Particle::Bank* bank_array, int64_t* bank_size, int64_t bank_capacity);
+  std::vector<Particle::Bank>& bank);
 
 int sample_element(Particle* p);
 
diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h
index cdf1f7089c..f246cb6ecb 100644
--- a/include/openmc/physics_mg.h
+++ b/include/openmc/physics_mg.h
@@ -8,6 +8,8 @@
 #include "openmc/particle.h"
 #include "openmc/nuclide.h"
 
+#include <vector>
+
 namespace openmc {
 
 //! \brief samples particle behavior after a collision event.
@@ -31,12 +33,9 @@ scatter(Particle* p);
 //! \brief Determines the average total, prompt and delayed neutrons produced
 //! from fission and creates the appropriate bank sites.
 //! \param p Particle to operate on
-//! \param bank_array The particle bank to populate
-//! \param size_bank Number of particles currently in the bank
-//! \param bank_array_size Allocated size of the bank
+//! \param bank The particle bank to populate
 void
-create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
-     int64_t bank_array_size);
+create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank);
 
 //! \brief Handles an absorption event
 //! \param p Particle to operate on
diff --git a/src/bank.cpp b/src/bank.cpp
index 8b7b436889..fc48a77615 100644
--- a/src/bank.cpp
+++ b/src/bank.cpp
@@ -16,8 +16,6 @@ namespace openmc {
 
 namespace simulation {
 
-int64_t n_bank;
-
 std::vector<Particle::Bank> source_bank;
 std::vector<Particle::Bank> fission_bank;
 #ifdef _OPENMP
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index c995a02d7d..1701955295 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -23,6 +23,7 @@
 #include <algorithm> // for min
 #include <array>
 #include <cmath> // for sqrt, abs, pow
+#include <iterator> // for back_inserter
 #include <string>
 
 namespace openmc {
@@ -81,21 +82,22 @@ void synchronize_bank()
 
 #ifdef OPENMC_MPI
   int64_t start = 0;
-  MPI_Exscan(&simulation::n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
+  int64_t n_bank = simulation::fission_bank.size();
+  MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
 
   // While we would expect the value of start on rank 0 to be 0, the MPI
   // standard says that the receive buffer on rank 0 is undefined and not
   // significant
   if (mpi::rank == 0) start = 0;
 
-  int64_t finish = start + simulation::n_bank;
+  int64_t finish = start + simulation::fission_bank.size();
   int64_t total = finish;
   MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
 
 #else
   int64_t start  = 0;
-  int64_t finish = simulation::n_bank;
-  int64_t total  = simulation::n_bank;
+  int64_t finish = simulation::fission_bank.size();
+  int64_t total  = simulation::fission_bank.size();
 #endif
 
   // If there are not that many particles per generation, it's possible that no
@@ -103,7 +105,7 @@ void synchronize_bank()
   // extra logic to treat this circumstance, we really want to ensure the user
   // runs enough particles to avoid this in the first place.
 
-  if (simulation::n_bank == 0) {
+  if (simulation::fission_bank.empty()) {
     fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
   }
 
@@ -134,21 +136,21 @@ void synchronize_bank()
   int64_t index_temp = 0;
   std::vector<Particle::Bank> temp_sites(3*simulation::work);
 
-  for (int64_t i = 0; i < simulation::n_bank; ++i) {
+  for (const auto& site : simulation::fission_bank) {
     // If there are less than n_particles particles banked, automatically add
     // int(n_particles/total) sites to temp_sites. For example, if you need
     // 1000 and 300 were banked, this would add 3 source sites per banked site
     // and the remaining 100 would be randomly sampled.
     if (total < settings::n_particles) {
       for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
-        temp_sites[index_temp] = simulation::fission_bank[i];
+        temp_sites[index_temp] = site;
         ++index_temp;
       }
     }
 
     // Randomly sample sites needed
     if (prn() < p_sample) {
-      temp_sites[index_temp] = simulation::fission_bank[i];
+      temp_sites[index_temp] = site;
       ++index_temp;
     }
   }
@@ -189,7 +191,7 @@ void synchronize_bank()
       // fission bank
       sites_needed = settings::n_particles - finish;
       for (int i = 0; i < sites_needed; ++i) {
-        int i_bank = simulation::n_bank - sites_needed + i;
+        int i_bank = simulation::fission_bank.size() - sites_needed + i;
         temp_sites[index_temp] = simulation::fission_bank[i_bank];
         ++index_temp;
       }
@@ -346,38 +348,30 @@ void calculate_average_keff()
 #ifdef _OPENMP
 void join_bank_from_threads()
 {
-  // Initialize the total number of fission bank sites
-  int64_t total = 0;
-
-#pragma omp parallel
+  #pragma omp parallel
   {
     // Copy thread fission bank sites to one shared copy
-#pragma omp for ordered schedule(static)
+    #pragma omp for ordered schedule(static)
     for (int i = 0; i < simulation::n_threads; ++i) {
-#pragma omp ordered
+      #pragma omp ordered
       {
         std::copy(
-          &simulation::fission_bank[0],
-          &simulation::fission_bank[0] + simulation::n_bank,
-          &simulation::master_fission_bank[total]
+          simulation::fission_bank.cbegin(),
+          simulation::fission_bank.cend(),
+          std::back_inserter(simulation::master_fission_bank)
         );
-        total += simulation::n_bank;
       }
     }
 
     // Make sure all threads have made it to this point
-#pragma omp barrier
+    #pragma omp barrier
 
     // Now copy the shared fission bank sites back to the master thread's copy.
     if (simulation::thread_id == 0) {
-      simulation::n_bank = total;
-      std::copy(
-        &simulation::master_fission_bank[0],
-        &simulation::master_fission_bank[0] + simulation::n_bank,
-        &simulation::fission_bank[0]
-      );
+      simulation::fission_bank = simulation::master_fission_bank;
+      simulation::master_fission_bank.clear();
     } else {
-      simulation::n_bank = 0;
+      simulation::fission_bank.clear();
     }
   }
 }
@@ -537,8 +531,8 @@ void shannon_entropy()
 
   // Get source weight in each mesh bin
   bool sites_outside;
-  xt::xtensor<double, 1> p = m->count_sites(simulation::n_bank,
-    simulation::fission_bank.data(), 0, nullptr, &sites_outside);
+  xt::xtensor<double, 1> p = m->count_sites(simulation::fission_bank,
+    &sites_outside);
 
   // display warning message if there were sites outside entropy box
   if (sites_outside) {
@@ -577,8 +571,8 @@ void ufs_count_sites()
   } else {
     // count number of source sites in each ufs mesh cell
     bool sites_outside;
-    simulation::source_frac = m->count_sites(simulation::work,
-      simulation::source_bank.data(), 0, nullptr, &sites_outside);
+    simulation::source_frac = m->count_sites(simulation::source_bank,
+      &sites_outside);
 
     // Check for sites outside of the mesh
     if (mpi::master && sites_outside) {
diff --git a/src/mesh.cpp b/src/mesh.cpp
index b75a0aad1f..ddfa272a6d 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -732,25 +732,21 @@ void RegularMesh::to_hdf5(hid_t group) const
   close_group(mesh_group);
 }
 
-xt::xarray<double> RegularMesh::count_sites(int64_t n, const Particle::Bank* bank,
-  int n_energy, const double* energies, bool* outside) const
+xt::xarray<double>
+RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
+  bool* outside) const
 {
   // Determine shape of array for counts
   std::size_t m = xt::prod(shape_)();
-  std::vector<std::size_t> shape;
-  if (n_energy > 0) {
-    shape = {m, static_cast<std::size_t>(n_energy - 1)};
-  } else {
-    shape = {m};
-  }
+  std::vector<std::size_t> shape = {m};
 
   // Create array of zeros
   xt::xarray<double> cnt {shape, 0.0};
   bool outside_ = false;
 
-  for (int64_t i = 0; i < n; ++i) {
+  for (const auto& site : bank) {
     // determine scoring bin for entropy mesh
-    int mesh_bin = get_bin(bank[i].r);
+    int mesh_bin = get_bin(site.r);
 
     // if outside mesh, skip particle
     if (mesh_bin < 0) {
@@ -758,19 +754,8 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Particle::Bank* ban
       continue;
     }
 
-    if (n_energy > 0) {
-      double E = bank[i].E;
-      if (E >= energies[0] && E <= energies[n_energy - 1]) {
-        // determine energy bin
-        int e_bin = lower_bound_index(energies, energies + n_energy, E);
-
-        // Add to appropriate bin
-        cnt(mesh_bin, e_bin) += bank[i].wgt;
-      }
-    } else {
-      // Add to appropriate bin
-      cnt(mesh_bin) += bank[i].wgt;
-    }
+    // Add to appropriate bin
+    cnt(mesh_bin) += site.wgt;
   }
 
   // Create copy of count data
diff --git a/src/particle.cpp b/src/particle.cpp
index c6f73f4f12..92045e2cdb 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -75,17 +75,14 @@ Particle::clear()
 void
 Particle::create_secondary(Direction u, double E, Type type)
 {
-  if (n_secondary_ == MAX_SECONDARY) {
-    fatal_error("Too many secondary particles created.");
-  }
+  secondary_bank_.emplace_back();
 
-  int64_t n = n_secondary_;
-  secondary_bank_[n].particle = type;
-  secondary_bank_[n].wgt = wgt_;
-  secondary_bank_[n].r = this->r();
-  secondary_bank_[n].u = u;
-  secondary_bank_[n].E = settings::run_CE ? E : g_;
-  ++n_secondary_;
+  auto& bank {secondary_bank_.back()};
+  bank.particle = type;
+  bank.wgt = wgt_;
+  bank.r = this->r();
+  bank.u = u;
+  bank.E = settings::run_CE ? E : g_;
 }
 
 void
@@ -360,10 +357,10 @@ Particle::transport()
     // Check for secondary particles if this particle is dead
     if (!alive_) {
       // If no secondary particles, break out of event loop
-      if (n_secondary_ == 0) break;
+      if (secondary_bank_.empty()) break;
 
-      this->from_source(&secondary_bank_[n_secondary_ - 1]);
-      --n_secondary_;
+      this->from_source(&secondary_bank_.back());
+      secondary_bank_.pop_back();
       n_event = 0;
 
       // Enter new particle in particle track file
diff --git a/src/physics.cpp b/src/physics.cpp
index 4cc9b12a03..7d44aa17c6 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -91,12 +91,10 @@ void sample_neutron_reaction(Particle* p)
   if (nuc->fissionable_) {
     Reaction* rx = sample_fission(i_nuclide, p);
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(),
-        &simulation::n_bank, simulation::fission_bank.size());
+      create_fission_sites(p, i_nuclide, rx, simulation::fission_bank);
     } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE &&
       settings::create_fission_neutrons) {
-      create_fission_sites(p, i_nuclide, rx, p->secondary_bank_,
-        &p->n_secondary_, MAX_SECONDARY);
+      create_fission_sites(p, i_nuclide, rx, p->secondary_bank_);
     }
   }
 
@@ -137,10 +135,8 @@ void sample_neutron_reaction(Particle* p)
 
 void
 create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
-  Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_capacity)
+  std::vector<Particle::Bank>& bank)
 {
-  // TODO: Heat generation from fission
-
   // If uniform fission source weighting is turned on, we increase or decrease
   // the expected number of fission sites produced
   double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
@@ -153,47 +149,30 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
   int nu = static_cast<int>(nu_t);
   if (prn() <= (nu_t - nu)) ++nu;
 
-  // Check for the bank size getting hit. For fixed source calculations, this
-  // is a fatal error; for eigenvalue calculations, it just means that k-eff
-  // was too high for a single batch.
-  if (*size_bank + nu > bank_capacity) {
-    if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
-      throw std::runtime_error{"Secondary particle bank size limit reached."
-           " If you are running a subcritical multiplication problem,"
-           " k-effective may be too close to one."};
-    } else {
-      if (mpi::master) {
-        warning("Maximum number of sites in fission bank reached. This can"
-          " result in irreproducible results using different numbers of"
-          " processes/threads.");
-      }
-    }
-  }
-
   // Begin banking the source neutrons
   // First, if our bank is full then don't continue
-  if (nu == 0 || *size_bank == bank_capacity) return;
+  if (nu == 0) return;
 
   // Initialize the counter of delayed neutrons encountered for each delayed
   // group.
   double nu_d[MAX_DELAYED_GROUPS] = {0.};
 
   p->fission_ = true;
-  for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) {
+  for (int i = 0; i < nu; ++i) {
+    // Create new bank site and get reference to last element
+    bank.emplace_back();
+    auto& site {bank.back()};
+
     // Bank source neutrons by copying the particle data
-    bank_array[i].r = p->r();
-
-    // Set that the bank particle is a neutron
-    bank_array[i].particle = Particle::Type::neutron;
-
-    // Set the weight of the fission bank site
-    bank_array[i].wgt = 1. / weight;
+    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_, &bank_array[i]);
+    sample_fission_neutron(i_nuclide, rx, p->E_, &site);
 
     // Set the delayed group on the particle as well
-    p->delayed_group_ = bank_array[i].delayed_group;
+    p->delayed_group_ = site.delayed_group;
 
     // Increment the number of neutrons born delayed
     if (p->delayed_group_ > 0) {
@@ -201,9 +180,6 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
     }
   }
 
-  // Increment number of bank sites
-  *size_bank = std::min(*size_bank + nu, bank_capacity);
-
   // Store the total weight banked for analog fission tallies
   p->n_bank_ = nu;
   p->wgt_bank_ = nu / weight;
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index affe0ec197..2a6b700f03 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -48,13 +48,10 @@ sample_reaction(Particle* p)
 
   if (model::materials[p->material_]->fissionable_) {
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      create_fission_sites(
-        p, simulation::fission_bank.data(), &simulation::n_bank,
-        simulation::fission_bank.size());
+      create_fission_sites(p, simulation::fission_bank);
     } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
                (settings::create_fission_neutrons)) {
-      create_fission_sites(p, p->secondary_bank_, &(p->n_secondary_),
-                           MAX_SECONDARY);
+      create_fission_sites(p, p->secondary_bank_);
     }
   }
 
@@ -102,11 +99,8 @@ scatter(Particle* p)
 }
 
 void
-create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
-     int64_t bank_array_size)
+create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
 {
-  // TODO: Heat generation from fission
-
   // If uniform fission source weighting is turned on, we increase or decrease
   // the expected number of fission sites produced
   double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
@@ -121,62 +115,42 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank
     nu++;
   }
 
-  // Check for the bank size getting hit. For fixed source calculations, this
-  // is a fatal error; for eigenvalue calculations, it just means that k-eff
-  // was too high for a single batch.
-  if (*size_bank + nu > bank_array_size) {
-    if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
-      throw std::runtime_error{"Secondary particle bank size limit reached."
-           " If you are running a subcritical multiplication problem,"
-           " k-effective may be too close to one."};
-    } else {
-      if (mpi::master) {
-        std::stringstream msg;
-        msg << "Maximum number of sites in fission bank reached. This can"
-             " result in irreproducible results using different numbers of"
-             " processes/threads.";
-        warning(msg);
-      }
-    }
-  }
-
   // Begin banking the source neutrons
   // First, if our bank is full then don't continue
-  if ((nu == 0) || (*size_bank == bank_array_size)) return;
+  if (nu == 0) return;
 
   // Initialize the counter of delayed neutrons encountered for each delayed
   // group.
   double nu_d[MAX_DELAYED_GROUPS] = {0.};
 
   p->fission_ = true;
-  for (size_t i = static_cast<size_t>(*size_bank);
-       i < static_cast<size_t>(std::min(*size_bank + nu, bank_array_size)); i++) {
+  for (int i = 0; i < nu; ++i) {
+    // Create new bank site and get reference to last element
+    bank.emplace_back();
+    auto& site {bank.back()};
+
     // Bank source neutrons by copying the particle data
-    bank_array[i].r = p->r();
-
-    // Set that the bank particle is a neutron
-    bank_array[i].particle = Particle::Type::neutron;
-
-    // Set the weight of the fission bank site
-    bank_array[i].wgt = 1. / weight;
+    site.r = p->r();
+    site.particle = Particle::Type::neutron;
+    site.wgt = 1. / weight;
 
     // Sample the cosine of the angle, assuming fission neutrons are emitted
     // isotropically
-    double mu = 2. * prn() - 1.;
+    double mu = 2.*prn() - 1.;
 
     // Sample the azimuthal angle uniformly in [0, 2.pi)
     double phi = 2. * PI * prn();
-    bank_array[i].u.x = mu;
-    bank_array[i].u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
-    bank_array[i].u.z = std::sqrt(1. - mu * mu) * std::sin(phi);
+    site.u.x = mu;
+    site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
+    site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi);
 
     // Sample secondary energy distribution for the fission reaction and set
     // the energy in the fission bank
     int dg;
     int gout;
     data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout);
-    bank_array[i].E = gout + 1;
-    bank_array[i].delayed_group = dg + 1;
+    site.E = gout + 1;
+    site.delayed_group = dg + 1;
 
     // Set the delayed group on the particle as well
     p->delayed_group_ = dg + 1;
@@ -187,9 +161,6 @@ create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank
     }
   }
 
-  // Increment number of bank sites
-  *size_bank = std::min(*size_bank + nu, bank_array_size);
-
   // Store the total weight banked for analog fission tallies
   p->n_bank_ = nu;
   p->wgt_bank_ = nu / weight;
diff --git a/src/simulation.cpp b/src/simulation.cpp
index b4db5e0e6d..fca0c85282 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -279,14 +279,14 @@ void allocate_banks()
     {
       simulation::thread_id = omp_get_thread_num();
       if (simulation::thread_id == 0) {
-        simulation::fission_bank.resize(3*simulation::work);
+        simulation::fission_bank.reserve(3*simulation::work);
       } else {
-        simulation::fission_bank.resize(3*simulation::work / simulation::n_threads);
+        simulation::fission_bank.reserve(3*simulation::work / simulation::n_threads);
       }
     }
-    simulation::master_fission_bank.resize(3*simulation::work);
+    simulation::master_fission_bank.reserve(3*simulation::work);
 #else
-    simulation::fission_bank.resize(3*simulation::work);
+    simulation::fission_bank.reserve(3*simulation::work);
 #endif
   }
 }
@@ -386,8 +386,8 @@ void finalize_batch()
 void initialize_generation()
 {
   if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-    // Reset number of fission bank sites
-    simulation::n_bank = 0;
+    // Clear out the fission bank
+    simulation::fission_bank.clear();
 
     // Count source sites if using uniform fission source weighting
     if (settings::ufs_on) ufs_count_sites();
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index de851ea3cb..7e235fa458 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -192,7 +192,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
 
   // loop over number of particles banked
   for (auto i = 0; i < p->n_bank_; ++i) {
-    auto i_bank = simulation::n_bank - p->n_bank_ + i;
+    auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
     const auto& bank = simulation::fission_bank[i_bank];
 
     // get the delayed group
@@ -813,7 +813,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
           // contribution to the fission bank to the score.
           score = 0.;
           for (auto i = 0; i < p->n_bank_; ++i) {
-            auto i_bank = simulation::n_bank - p->n_bank_ + i;
+            auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
             const auto& bank = simulation::fission_bank[i_bank];
             auto g = bank.delayed_group;
             if (g != 0) {
@@ -1785,7 +1785,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           // contribution to the fission bank to the score.
           score = 0.;
           for (auto i = 0; i < p->n_bank_; ++i) {
-            auto i_bank = simulation::n_bank - p->n_bank_ + i;
+            auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i;
             const auto& bank = simulation::fission_bank[i_bank];
             auto g = bank.delayed_group;
             if (g != 0) {