diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 6ac8aabae..00476e3ef 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -218,6 +218,10 @@ public:
   //! create a particle restart HDF5 file
   void write_restart() const;
 
+  //! Gets the pointer to the particle's current PRN seed
+  uint64_t* current_seed();
+  const uint64_t* current_seed() const;
+
   //==========================================================================
   // Data members
 
diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h
index 350fa7d24..5521f1df0 100644
--- a/include/openmc/random_lcg.h
+++ b/include/openmc/random_lcg.h
@@ -43,11 +43,13 @@ double future_prn(int64_t n, uint64_t prn_seed);
 //==============================================================================
 //! Set a RNG seed to a unique value based on a unique particle ID by striding
 //! the seed.
-//! @param prn_seeds Pseudorandom number seed array
 //! @param id The particle ID
+//! @param offset The offset from the master seed to be used (e.g., for creating
+//! different streams)
+//! @return The initialized seed value
 //==============================================================================
 
-void init_seed(int64_t id, uint64_t* prn_seeds, int offset );
+uint64_t init_seed(int64_t id, int offset);
 
 //==============================================================================
 //! Set the RNG seeds to unique values based on the ID of the particle. This
@@ -57,7 +59,7 @@ void init_seed(int64_t id, uint64_t* prn_seeds, int offset );
 //! @param id The particle ID
 //==============================================================================
 
-void init_particle_seeds(int64_t id, uint64_t* prn_seeds );
+void init_particle_seeds(int64_t id, uint64_t* prn_seeds);
 
 //==============================================================================
 //! Advance the random number seed 'n' times from the current seed. This
diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp
index 85733eccc..e7476124d 100644
--- a/src/bremsstrahlung.cpp
+++ b/src/bremsstrahlung.cpp
@@ -65,7 +65,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
   double y = std::exp(y_l + (y_r - y_l)*f);
 
   // Sample number of secondary bremsstrahlung photons
-  int n = y + prn(p.prn_seeds_ + p.stream_);
+  int n = y + prn(p.current_seed());
 
   *E_lost = 0.0;
   if (n == 0) return;
@@ -73,7 +73,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
   // Sample index of the tabulated PDF in the energy grid, j or j+1
   double c_max;
   int i_e;
-  if (prn(p.prn_seeds_ + p.stream_) <= f || j == 0) {
+  if (prn(p.current_seed()) <= f || j == 0) {
     i_e = j + 1;
 
     // Interpolate the maximum value of the CDF at the incoming particle
@@ -94,7 +94,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
   for (int i = 0; i < n; ++i) {
     // Generate a random number r and determine the index i for which
     // cdf(i) <= r*cdf,max <= cdf(i+1)
-    double c = prn(p.prn_seeds_ + p.stream_)*c_max;
+    double c = prn(p.current_seed())*c_max;
     int i_w = lower_bound_index(&mat->cdf(i_e, 0), &mat->cdf(i_e, 0) + i_e, c);
 
     // Sample the photon energy
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index cd0699a32..110f0100e 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -114,16 +114,13 @@ void synchronize_bank()
     fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
   }
 
-  // Create pseudorandom number seed
-  uint64_t prn_seed;
-
   // Make sure all processors start at the same point for random sampling. Then
   // skip ahead in the sequence using the starting index in the 'global'
   // fission bank for each processor.
 
   int64_t id = simulation::total_gen + overall_generation();
-  init_seed(id, &prn_seed, STREAM_TRACKING);
-  advance_prn_seed(start, &prn_seed);
+  uint64_t seed = init_seed(id, STREAM_TRACKING);
+  advance_prn_seed(start, &seed);
 
   // Determine how many fission sites we need to sample from the source bank
   // and the probability for selecting a site.
@@ -158,7 +155,7 @@ void synchronize_bank()
     }
 
     // Randomly sample sites needed
-    if (prn(&prn_seed) < p_sample) {
+    if (prn(&seed) < p_sample) {
       temp_sites[index_temp] = site;
       ++index_temp;
     }
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index bc0fed2e0..86f2ee140 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -581,7 +581,7 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle
 
       // Randomly sample between temperature i and i+1
       f = (kT - kTs_[i_temp]) / (kTs_[i_temp + 1] - kTs_[i_temp]);
-      if (f > prn(p.prn_seeds_ + p.stream_)) ++i_temp;
+      if (f > prn(p.current_seed())) ++i_temp;
       break;
     }
 
@@ -720,7 +720,7 @@ void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p)
   int i_temp;
   double elastic;
   double inelastic;
-  data::thermal_scatt[i_sab]->calculate_xs(p.E_, p.sqrtkT_, &i_temp, &elastic, &inelastic, p.prn_seeds_ + p.stream_);
+  data::thermal_scatt[i_sab]->calculate_xs(p.E_, p.sqrtkT_, &i_temp, &elastic, &inelastic, p.current_seed());
 
   // Store the S(a,b) cross sections.
   micro.thermal = sab_frac * (elastic + inelastic);
@@ -759,7 +759,7 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
   p.stream_ = STREAM_URR_PTABLE;
   //TODO: to maintain the same random number stream as the Fortran code this
   //replaces, the seed is set with i_nuclide_ + 1 instead of i_nuclide_
-  double r = future_prn(static_cast<int64_t>(i_nuclide_ + 1), p.prn_seeds_[p.stream_]);
+  double r = future_prn(static_cast<int64_t>(i_nuclide_ + 1), *p.current_seed());
   p.stream_ = STREAM_TRACKING;
 
   int i_low = 0;
diff --git a/src/particle.cpp b/src/particle.cpp
index d3db2e84b..6fb41f295 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -228,7 +228,7 @@ Particle::transport()
     } else if (macro_xs_.total == 0.0) {
       d_collision = INFINITY;
     } else {
-      d_collision = -std::log(prn(prn_seeds_ + stream_)) / macro_xs_.total;
+      d_collision = -std::log(prn(this->current_seed())) / macro_xs_.total;
     }
 
     // Select smaller of the two distances
@@ -461,7 +461,7 @@ Particle::cross_surface()
 
     Direction u = (surf->bc_ == BC_REFLECT) ?
       surf->reflect(this->r(), this->u()) :
-      surf->diffuse_reflect(this->r(), this->u(), prn_seeds_ + stream_);
+      surf->diffuse_reflect(this->r(), this->u(), this->current_seed());
 
     // Make sure new particle direction is normalized
     this->u() = u / u.norm();
@@ -678,4 +678,7 @@ Particle::write_restart() const
   } // #pragma omp critical
 }
 
+uint64_t* Particle::current_seed() {return prn_seeds_ + stream_;}
+const uint64_t* Particle::current_seed() const {return prn_seeds_ + stream_;}
+
 } // namespace openmc
diff --git a/src/photon.cpp b/src/photon.cpp
index 6cf43cbc5..066e182fc 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -649,8 +649,8 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
 {
   // If no transitions, assume fluorescent photon from captured free electron
   if (shell.n_transitions == 0) {
-    double mu = 2.0*prn(p.prn_seeds_ + p.stream_) - 1.0;
-    double phi = 2.0*PI*prn(p.prn_seeds_ + p.stream_);
+    double mu = 2.0*prn(p.current_seed()) - 1.0;
+    double phi = 2.0*PI*prn(p.current_seed());
     Direction u;
     u.x = mu;
     u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
@@ -661,7 +661,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
   }
 
   // Sample transition
-  double rn = prn(p.prn_seeds_ + p.stream_);
+  double rn = prn(p.current_seed());
   double c = 0.0;
   int i_transition;
   for (i_transition = 0; i_transition < shell.n_transitions; ++i_transition) {
@@ -674,8 +674,8 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
   int secondary = shell.transition_subshells(i_transition, 1);
 
   // Sample angle isotropically
-  double mu = 2.0*prn(p.prn_seeds_ + p.stream_) - 1.0;
-  double phi = 2.0*PI*prn(p.prn_seeds_ + p.stream_);
+  double mu = 2.0*prn(p.current_seed()) - 1.0;
+  double phi = 2.0*PI*prn(p.current_seed());
   Direction u;
   u.x = mu;
   u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
diff --git a/src/physics.cpp b/src/physics.cpp
index 53e4a5ac0..7b86259fe 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -134,7 +134,7 @@ void sample_neutron_reaction(Particle* p)
   // Advance URR seed stream 'N' times after energy changes
   if (p->E_ != p->E_last_) {
     p->stream_ = STREAM_URR_PTABLE;
-    advance_prn_seed(data::nuclides.size(), p->prn_seeds_ + p->stream_);
+    advance_prn_seed(data::nuclides.size(), p->current_seed());
     p->stream_ = STREAM_TRACKING;
   }
 
@@ -159,7 +159,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 
   // Sample the number of neutrons produced
   int nu = static_cast<int>(nu_t);
-  if (prn(p->prn_seeds_ + p->stream_) <= (nu_t - nu)) ++nu;
+  if (prn(p->current_seed()) <= (nu_t - nu)) ++nu;
 
   // Begin banking the source neutrons
   // First, if our bank is full then don't continue
@@ -181,7 +181,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
     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->current_seed());
 
     // Set the delayed group on the particle as well
     p->delayed_group_ = site.delayed_group;
@@ -223,13 +223,13 @@ void sample_photon_reaction(Particle* p)
   // For tallying purposes, this routine might be called directly. In that
   // case, we need to sample a reaction via the cutoff variable
   double prob = 0.0;
-  double cutoff = prn(p->prn_seeds_ + p->stream_) * micro.total;
+  double cutoff = prn(p->current_seed()) * micro.total;
 
   // Coherent (Rayleigh) scattering
   prob += micro.coherent;
   if (prob > cutoff) {
-    double mu = element.rayleigh_scatter(alpha, p->prn_seeds_ + p->stream_);
-    p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_ + p->stream_);
+    double mu = element.rayleigh_scatter(alpha, p->current_seed());
+    p->u() = rotate_angle(p->u(), mu, nullptr, p->current_seed());
     p->event_ = EVENT_SCATTER;
     p->event_mt_ = COHERENT;
     return;
@@ -240,7 +240,7 @@ void sample_photon_reaction(Particle* p)
   if (prob > cutoff) {
     double alpha_out, mu;
     int i_shell;
-    element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell, p->prn_seeds_ + p->stream_);
+    element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell, p->current_seed());
 
     // Determine binding energy of shell. The binding energy is 0.0 if
     // doppler broadening is not used.
@@ -252,13 +252,13 @@ void sample_photon_reaction(Particle* p)
     }
 
     // Create Compton electron
-    double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
+    double phi = 2.0*PI*prn(p->current_seed());
     double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
     int electron = static_cast<int>(Particle::Type::electron);
     if (E_electron >= settings::energy_cutoff[electron]) {
       double mu_electron = (alpha - alpha_out*mu)
         / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
-      Direction u = rotate_angle(p->u(), mu_electron, &phi, p->prn_seeds_ + p->stream_);
+      Direction u = rotate_angle(p->u(), mu_electron, &phi, p->current_seed());
       p->create_secondary(u, E_electron, Particle::Type::electron);
     }
 
@@ -272,7 +272,7 @@ void sample_photon_reaction(Particle* p)
 
     phi += PI;
     p->E_ = alpha_out*MASS_ELECTRON_EV;
-    p->u() = rotate_angle(p->u(), mu, &phi, p->prn_seeds_ + p->stream_);
+    p->u() = rotate_angle(p->u(), mu, &phi, p->current_seed());
     p->event_ = EVENT_SCATTER;
     p->event_mt_ = INCOHERENT;
     return;
@@ -304,8 +304,8 @@ void sample_photon_reaction(Particle* p)
         // model in Serpent 2" by Toni Kaltiaisenaho
         double mu;
         while (true) {
-          double r = prn(p->prn_seeds_ + p->stream_);
-          if (4.0*(1.0 - r)*r >= prn(p->prn_seeds_ + p->stream_)) {
+          double r = prn(p->current_seed());
+          if (4.0*(1.0 - r)*r >= prn(p->current_seed())) {
             double rel_vel = std::sqrt(E_electron * (E_electron +
               2.0*MASS_ELECTRON_EV)) / (E_electron + MASS_ELECTRON_EV);
             mu = (2.0*r + rel_vel - 1.0) / (2.0*rel_vel*r - rel_vel + 1.0);
@@ -313,7 +313,7 @@ void sample_photon_reaction(Particle* p)
           }
         }
 
-        double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
+        double phi = 2.0*PI*prn(p->current_seed());
         Direction u;
         u.x = mu;
         u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
@@ -341,14 +341,14 @@ void sample_photon_reaction(Particle* p)
     double E_electron, E_positron;
     double mu_electron, mu_positron;
     element.pair_production(alpha, &E_electron, &E_positron,
-      &mu_electron, &mu_positron, p->prn_seeds_ + p->stream_);
+      &mu_electron, &mu_positron, p->current_seed());
 
     // Create secondary electron
-    Direction u = rotate_angle(p->u(), mu_electron, nullptr, p->prn_seeds_ + p->stream_);
+    Direction u = rotate_angle(p->u(), mu_electron, nullptr, p->current_seed());
     p->create_secondary(u, E_electron, Particle::Type::electron);
 
     // Create secondary positron
-    u = rotate_angle(p->u(), mu_positron, nullptr, p->prn_seeds_ + p->stream_);
+    u = rotate_angle(p->u(), mu_positron, nullptr, p->current_seed());
     p->create_secondary(u, E_positron, Particle::Type::positron);
 
     p->event_ = EVENT_ABSORB;
@@ -382,8 +382,8 @@ void sample_positron_reaction(Particle* p)
   }
 
   // Sample angle isotropically
-  double mu = 2.0*prn(p->prn_seeds_ + p->stream_) - 1.0;
-  double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
+  double mu = 2.0*prn(p->current_seed()) - 1.0;
+  double phi = 2.0*PI*prn(p->current_seed());
   Direction u;
   u.x = mu;
   u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
@@ -401,7 +401,7 @@ void sample_positron_reaction(Particle* p)
 int sample_nuclide(Particle* p)
 {
   // Sample cumulative distribution function
-  double cutoff = prn(p->prn_seeds_ + p->stream_) * p->macro_xs_.total;
+  double cutoff = prn(p->current_seed()) * p->macro_xs_.total;
 
   // Get pointers to nuclide/density arrays
   const auto& mat {model::materials[p->material_]};
@@ -426,7 +426,7 @@ int sample_nuclide(Particle* p)
 int sample_element(Particle* p)
 {
   // Sample cumulative distribution function
-  double cutoff = prn(p->prn_seeds_ + p->stream_) * p->macro_xs_.total;
+  double cutoff = prn(p->current_seed()) * p->macro_xs_.total;
 
   // Get pointers to elements, densities
   const auto& mat {model::materials[p->material_]};
@@ -479,7 +479,7 @@ Reaction* sample_fission(int i_nuclide, Particle* p)
   int i_temp = p->neutron_xs_[i_nuclide].index_temp;
   int i_grid = p->neutron_xs_[i_nuclide].index_grid;
   double f = p->neutron_xs_[i_nuclide].interp_factor;
-  double cutoff = prn(p->prn_seeds_ + p->stream_) * p->neutron_xs_[i_nuclide].fission;
+  double cutoff = prn(p->current_seed()) * p->neutron_xs_[i_nuclide].fission;
   double prob = 0.0;
 
   // Loop through each partial fission reaction type
@@ -506,7 +506,7 @@ void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product
   int i_temp = p->neutron_xs_[i_nuclide].index_temp;
   int i_grid = p->neutron_xs_[i_nuclide].index_grid;
   double f = p->neutron_xs_[i_nuclide].interp_factor;
-  double cutoff = prn(p->prn_seeds_ + p->stream_) * p->neutron_xs_[i_nuclide].photon_prod;
+  double cutoff = prn(p->current_seed()) * p->neutron_xs_[i_nuclide].photon_prod;
   double prob = 0.0;
 
   // Loop through each reaction type
@@ -554,7 +554,7 @@ void absorption(Particle* p, int i_nuclide)
   } else {
     // See if disappearance reaction happens
     if (p->neutron_xs_[i_nuclide].absorption >
-        prn(p->prn_seeds_ + p->stream_) * p->neutron_xs_[i_nuclide].total) {
+        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_[
@@ -582,7 +582,7 @@ void scatter(Particle* p, int i_nuclide)
 
   // For tallying purposes, this routine might be called directly. In that
   // case, we need to sample a reaction via the cutoff variable
-  double cutoff = prn(p->prn_seeds_ + p->stream_) * (micro.total - micro.absorption);
+  double cutoff = prn(p->current_seed()) * (micro.total - micro.absorption);
   bool sampled = false;
 
   // Calculate elastic cross section if it wasn't precalculated
@@ -656,8 +656,8 @@ void scatter(Particle* p, int i_nuclide)
     int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
     if (mat->p0_[i_nuc_mat]) {
       // Sample isotropic-in-lab outgoing direction
-      double mu = 2.0*prn(p->prn_seeds_ + p->stream_) - 1.0;
-      double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
+      double mu = 2.0*prn(p->current_seed()) - 1.0;
+      double phi = 2.0*PI*prn(p->current_seed());
 
       // Change direction of particle
       p->u().x = mu;
@@ -684,7 +684,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   Direction v_t {};
   if (!p->neutron_xs_[i_nuclide].use_ptable) {
     v_t = sample_target_velocity(nuc.get(), p->E_, p->u(), v_n,
-      p->neutron_xs_[i_nuclide].elastic, kT, p->prn_seeds_ + p->stream_);
+      p->neutron_xs_[i_nuclide].elastic, kT, p->current_seed());
   }
 
   // Velocity of center-of-mass
@@ -702,9 +702,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   auto& d = rx.products_[0].distribution_[0];
   auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
   if (d_) {
-    mu_cm = d_->angle().sample(p->E_, p->prn_seeds_ + p->stream_);
+    mu_cm = d_->angle().sample(p->E_, p->current_seed());
   } else {
-    mu_cm = 2.0*prn(p->prn_seeds_ + p->stream_) - 1.0;
+    mu_cm = 2.0*prn(p->current_seed()) - 1.0;
   }
 
   // Determine direction cosines in CM
@@ -713,7 +713,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   // Rotate neutron velocity vector to new angle -- note that the speed of the
   // neutron in CM does not change in elastic scattering. However, the speed
   // will change when we convert back to LAB
-  v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p->prn_seeds_ + p->stream_);
+  v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p->current_seed());
 
   // Transform back to LAB frame
   v_n += v_cm;
@@ -742,11 +742,11 @@ void sab_scatter(int i_nuclide, int i_sab, Particle* p)
 
   // Sample energy and angle
   double E_out;
-  data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p->E_, &E_out, &p->mu_, p->prn_seeds_ + p->stream_);
+  data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p->E_, &E_out, &p->mu_, p->current_seed());
 
   // Set energy to outgoing, change direction of particle
   p->E_ = E_out;
-  p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->prn_seeds_ + p->stream_);
+  p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->current_seed());
 }
 
 Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
@@ -1050,7 +1050,7 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
   // sample outgoing energy and scattering cosine
   double E;
   double mu;
-  rx->products_[0].sample(E_in, E, mu, p->prn_seeds_ + p->stream_);
+  rx->products_[0].sample(E_in, E, mu, p->current_seed());
 
   // if scattering system is in center-of-mass, transfer cosine of scattering
   // angle and outgoing energy from CM to LAB
@@ -1076,7 +1076,7 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
   p->mu_ = mu;
 
   // change direction of particle
-  p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_ + p->stream_);
+  p->u() = rotate_angle(p->u(), mu, nullptr, p->current_seed());
 
   // evaluate yield
   double yield = (*rx->products_[0].yield_)(E_in);
@@ -1097,7 +1097,7 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
   double y_t = p->wgt_ * p->neutron_xs_[i_nuclide].photon_prod /
     p->neutron_xs_[i_nuclide].total;
   int y = static_cast<int>(y_t);
-  if (prn(p->prn_seeds_ + p->stream_) <= y_t - y) ++y;
+  if (prn(p->current_seed()) <= y_t - y) ++y;
 
   // Sample each secondary photon
   for (int i = 0; i < y; ++i) {
@@ -1110,10 +1110,10 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
     auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx];
     double E;
     double mu;
-    rx->products_[i_product].sample(p->E_, E, mu, p->prn_seeds_ + p->stream_);
+    rx->products_[i_product].sample(p->E_, E, mu, p->current_seed());
 
     // Sample the new direction
-    Direction u = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_ + p->stream_);
+    Direction u = rotate_angle(p->u(), mu, nullptr, p->current_seed());
 
     // Create the secondary photon
     p->create_secondary(u, E, Particle::Type::photon);
diff --git a/src/physics_common.cpp b/src/physics_common.cpp
index c26ebdfd2..1ea40a896 100644
--- a/src/physics_common.cpp
+++ b/src/physics_common.cpp
@@ -12,7 +12,7 @@ namespace openmc {
 void russian_roulette(Particle* p)
 {
   if (p->wgt_ < settings::weight_cutoff) {
-    if (prn(p->prn_seeds_ + p->stream_) < p->wgt_ / settings::weight_survive) {
+    if (prn(p->current_seed()) < p->wgt_ / settings::weight_survive) {
       p->wgt_ = settings::weight_survive;
       p->wgt_last_ = p->wgt_;
     } else {
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 3313f1034..f52e6c5bf 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -78,10 +78,10 @@ void
 scatter(Particle* p)
 {
   data::mg.macro_xs_[p->material_].sample_scatter(p->g_last_, p->g_, p->mu_,
-                                                  p->wgt_, p->prn_seeds_ + p->stream_);
+                                                  p->wgt_, p->current_seed());
 
   // Rotate the angle
-  p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->prn_seeds_ + p->stream_);
+  p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->current_seed());
 
   // Update energy value for downstream compatability (in tallying)
   p->E_ = data::mg.energy_bin_avg_[p->g_];
@@ -103,7 +103,7 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
 
   // Sample the number of neutrons produced
   int nu = static_cast<int>(nu_t);
-  if (prn(p->prn_seeds_ + p->stream_) <= (nu_t - int(nu_t))) {
+  if (prn(p->current_seed()) <= (nu_t - int(nu_t))) {
     nu++;
   }
 
@@ -128,10 +128,10 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
 
     // Sample the cosine of the angle, assuming fission neutrons are emitted
     // isotropically
-    double mu = 2.*prn(p->prn_seeds_ + p->stream_) - 1.;
+    double mu = 2.*prn(p->current_seed()) - 1.;
 
     // Sample the azimuthal angle uniformly in [0, 2.pi)
-    double phi = 2. * PI * prn(p->prn_seeds_ + p->stream_ );
+    double phi = 2. * PI * prn(p->current_seed() );
     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);
@@ -140,7 +140,7 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
     int dg;
     int gout;
     data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout,
-      p->prn_seeds_ + p->stream_);
+      p->current_seed());
     // Store the energy and delayed groups on the fission bank
     site.E = gout;
     // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest
@@ -180,7 +180,7 @@ absorption(Particle* p)
     global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission /
         p->macro_xs_.absorption;
   } else {
-    if (p->macro_xs_.absorption > prn(p->prn_seeds_ + p->stream_) * p->macro_xs_.total) {
+    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->macro_xs_.absorption;
diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp
index b5feb2ef0..dbae0f8c8 100644
--- a/src/random_lcg.cpp
+++ b/src/random_lcg.cpp
@@ -46,9 +46,9 @@ double future_prn(int64_t n, uint64_t prn_seed)
 // INIT_SEED
 //==============================================================================
 
-void init_seed(int64_t id, uint64_t* prn_seed, int offset)
+uint64_t init_seed(int64_t id, int offset)
 {
-  *prn_seed = future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + offset);
+  return future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + offset);
 }
 
 //==============================================================================
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
index f063e5685..d7c5ee63a 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -125,7 +125,7 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
     mu_out_(i, k+1) + f*(mu_out_(i+1, k+1) - mu_out_(i, k+1));
 
   // Smear cosine
-  mu += std::min(mu - mu_left, mu_right - mu)*(prn() - 0.5);
+  mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seed) - 0.5);
 
   // Energy doesn't change in elastic scattering
   E_out = E_in;
diff --git a/src/source.cpp b/src/source.cpp
index 483fd4b92..01d333649 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -270,11 +270,10 @@ void initialize_source()
       // initialize random number seed
       int64_t id = simulation::total_gen*settings::n_particles +
         simulation::work_index[mpi::rank] + i + 1;
-      uint64_t prn_seed;
-      init_seed(id, &prn_seed, STREAM_SOURCE);
+      uint64_t seed = init_seed(id, STREAM_SOURCE);
 
       // sample external source distribution
-      simulation::source_bank[i] = sample_external_source(&prn_seed);
+      simulation::source_bank[i] = sample_external_source(&seed);
     }
   }
 
@@ -331,11 +330,10 @@ void fill_source_bank_fixedsource()
       // initialize random number seed
       int64_t id = (simulation::total_gen + overall_generation()) *
         settings::n_particles + simulation::work_index[mpi::rank] + i + 1;
-      uint64_t prn_seed;
-      init_seed(id, &prn_seed, STREAM_SOURCE);
+      uint64_t seed = init_seed(id, STREAM_SOURCE);
 
       // sample external source distribution
-      simulation::source_bank[i] = sample_external_source(&prn_seed);
+      simulation::source_bank[i] = sample_external_source(&seed);
     }
   }
 }
diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp
index 7df5047f6..903a111a2 100644
--- a/src/volume_calc.cpp
+++ b/src/volume_calc.cpp
@@ -128,12 +128,11 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
       // Sample locations and count hits
       #pragma omp for
       for (size_t i = i_start; i < i_end; i++) {
-        uint64_t prn_seed;
         int64_t id = iterations * n_samples_ + i;
-        init_seed(id, &prn_seed, STREAM_VOLUME);
+        uint64_t seed = init_seed(id, STREAM_VOLUME);
 
         p.n_coord_ = 1;
-        Position xi {prn(&prn_seed), prn(&prn_seed), prn(&prn_seed)};
+        Position xi {prn(&seed), prn(&seed), prn(&seed)};
         p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
         p.u() = {0.5, 0.5, 0.5};