diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 1e06d9474..6ac8aabae 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -288,8 +288,8 @@ public:
   bool write_track_ {false};
 
   // Current PRNG state
-  uint64_t prn_seeds[N_STREAMS]; // current seeds
-  int      stream;               // current RNG stream
+  uint64_t prn_seeds_[N_STREAMS]; // current seeds
+  int      stream_;               // current RNG stream
 };
 
 } // namespace openmc
diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp
index 40cd28f7a..a99d118f0 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.prn_seeds_, p.stream_);
 
   *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.prn_seeds_, p.stream_) <= 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.prn_seeds_, p.stream_)*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/distribution_angle.cpp b/src/distribution_angle.cpp
index 41a316c11..59263afc8 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -62,7 +62,6 @@ AngleDistribution::AngleDistribution(hid_t group)
   }
 }
 
-//double AngleDistribution::sample(double E) const
 double AngleDistribution::sample(double E, uint64_t * prn_seeds, int stream) const
 {
   // Determine number of incoming energies
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index 9bf197378..fdfab2ef0 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -123,7 +123,7 @@ void synchronize_bank()
   // fission bank for each processor.
 
   set_particle_seed(simulation::total_gen + overall_generation(), prn_seeds);
-  advance_prn_seed(start, prn_seeds, stream); // TODO: What is starting stream here? Tracking? It doesn't appear to be enforced anywhere?
+  advance_prn_seed(start, prn_seeds, stream);
 
   // Determine how many fission sites we need to sample from the source bank
   // and the probability for selecting a site.
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index ccfd21189..b5da3f60a 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.prn_seeds_, p.stream_)) ++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.prn_seeds_, p.stream_);
 
   // Store the S(a,b) cross sections.
   micro.thermal = sab_frac * (elastic + inelastic);
@@ -756,11 +756,11 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
   // This guarantees the randomness and, at the same time, makes sure we
   // reuse random numbers for the same nuclide at different temperatures,
   // therefore preserving correlation of temperature in probability tables.
-  p.stream = STREAM_URR_PTABLE;
+  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);
-  p.stream = STREAM_TRACKING;
+  double r = future_prn(static_cast<int64_t>(i_nuclide_ + 1), p.prn_seeds_, p.stream_);
+  p.stream_ = STREAM_TRACKING;
 
   int i_low = 0;
   while (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;};
diff --git a/src/particle.cpp b/src/particle.cpp
index 10c58239f..98ee8920b 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -157,9 +157,9 @@ Particle::transport()
   while (true) {
     // Set the random number stream
     if (type_ == Particle::Type::neutron) {
-      stream = STREAM_TRACKING;
+      stream_ = STREAM_TRACKING;
     } else {
-      stream = STREAM_PHOTON;
+      stream_ = STREAM_PHOTON;
     }
 
     // Store pre-collision particle properties
@@ -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(prn_seeds_, stream_)) / 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(), prn_seeds_, stream_);
 
     // Make sure new particle direction is normalized
     this->u() = u / u.norm();
diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp
index 27ba1784d..17eb5a3e6 100644
--- a/src/particle_restart.cpp
+++ b/src/particle_restart.cpp
@@ -98,7 +98,7 @@ void run_particle_restart()
     throw std::runtime_error{"Unexpected run mode: " +
       std::to_string(previous_run_mode)};
   }
-  set_particle_seed(particle_seed, p.prn_seeds);
+  set_particle_seed(particle_seed, p.prn_seeds_);
 
   // Transport neutron
   p.transport();
diff --git a/src/photon.cpp b/src/photon.cpp
index cbb4133b2..803d89184 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.prn_seeds_, p.stream_) - 1.0;
+    double phi = 2.0*PI*prn(p.prn_seeds_, p.stream_);
     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.prn_seeds_, p.stream_);
   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.prn_seeds_, p.stream_) - 1.0;
+  double phi = 2.0*PI*prn(p.prn_seeds_, p.stream_);
   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 a17c0a8a4..28fc4341e 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -112,9 +112,9 @@ void sample_neutron_reaction(Particle* p)
 
   // Create secondary photons
   if (settings::photon_transport) {
-    p->stream = STREAM_PHOTON;
+    p->stream_ = STREAM_PHOTON;
     sample_secondary_photons(p, i_nuclide);
-	p->stream = STREAM_TRACKING;
+	p->stream_ = STREAM_TRACKING;
   }
 
   // If survival biasing is being used, the following subroutine adjusts the
@@ -133,9 +133,9 @@ 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);
-	p->stream = STREAM_TRACKING;
+	p->stream_ = STREAM_URR_PTABLE;
+    advance_prn_seed(data::nuclides.size(), p->prn_seeds_, p->stream_);
+	p->stream_ = STREAM_TRACKING;
   }
 
   // Play russian roulette if survival biasing is turned on
@@ -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->prn_seeds_, p->stream_) <= (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->prn_seeds_, p->stream_);
 
     // 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->prn_seeds_, p->stream_) * 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->prn_seeds_, p->stream_);
+    p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_, p->stream_);
     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->prn_seeds_, p->stream_);
 
     // 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->prn_seeds_, p->stream_);
     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->prn_seeds_, p->stream_);
       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->prn_seeds_, p->stream_);
     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->prn_seeds_, p->stream_);
+          if (4.0*(1.0 - r)*r >= prn(p->prn_seeds_, p->stream_)) {
             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->prn_seeds_, p->stream_);
         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->prn_seeds_, p->stream_);
 
     // 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->prn_seeds_, p->stream_);
     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->prn_seeds_, p->stream_);
     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->prn_seeds_, p->stream_) - 1.0;
+  double phi = 2.0*PI*prn(p->prn_seeds_, p->stream_);
   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->prn_seeds_, p->stream_) * 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->prn_seeds_, p->stream_) * 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->prn_seeds_, p->stream_) * 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->prn_seeds_, p->stream_) * 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->prn_seeds_, p->stream_) * 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->prn_seeds_, p->stream_) * (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->prn_seeds_, p->stream_) - 1.0;
+      double phi = 2.0*PI*prn(p->prn_seeds_, p->stream_);
 
       // 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->prn_seeds_, p->stream_);
   }
 
   // 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->prn_seeds_, p->stream_);
   } else {
-    mu_cm = 2.0*prn(p->prn_seeds, p->stream) - 1.0;
+    mu_cm = 2.0*prn(p->prn_seeds_, p->stream_) - 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->prn_seeds_, p->stream_);
 
   // 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->prn_seeds_, p->stream_);
 
   // 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->prn_seeds_, p->stream_);
 }
 
 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->prn_seeds_, p->stream_);
 
   // 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->prn_seeds_, p->stream_);
 
   // 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->prn_seeds_, p->stream_) <= 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->prn_seeds_, p->stream_);
 
     // 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->prn_seeds_, p->stream_);
 
     // 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 cfb4c01fe..e5b43e223 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->prn_seeds_, p->stream_) < 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 1d22a886b..6898e1e0e 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->prn_seeds_, p->stream_);
 
   // 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->prn_seeds_, p->stream_);
 
   // 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->prn_seeds_, p->stream_) <= (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->prn_seeds_, p->stream_) - 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->prn_seeds_, p->stream_ );
     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);
@@ -139,8 +139,8 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
     // Sample secondary energy distribution for the fission reaction
     int dg;
     int gout;
-    data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout, p->prn_seeds,
-        p->stream);
+    data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout, p->prn_seeds_,
+        p->stream_);
     // 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->prn_seeds_, p->stream_) * 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/secondary_thermal.cpp b/src/secondary_thermal.cpp
index 7c3af2993..fc4df50a6 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -32,7 +32,8 @@ CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs)
 { }
 
 void
-CoherentElasticAE::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
+CoherentElasticAE::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+  int stream) const
 {
   // Get index and interpolation factor for elastic grid
   int i;
@@ -65,7 +66,8 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group)
 }
 
 void
-IncoherentElasticAE::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
+IncoherentElasticAE::sample(double E_in, double& E_out, double& mu,
+  uint64_t * prn_seeds, int stream) const
 {
   // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
   double c = 2 * E_in * debye_waller_;
@@ -87,7 +89,8 @@ IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group,
 }
 
 void
-IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
+IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
+  uint64_t * prn_seeds, int stream) const
 {
   // Get index and interpolation factor for elastic grid
   int i;
@@ -125,7 +128,8 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group,
 }
 
 void
-IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
+IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
+  uint64_t * prn_seeds, int stream) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
@@ -228,7 +232,8 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
 }
 
 void
-IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
+IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
+  uint64_t * prn_seeds, int stream) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
index 1db62bf2f..f4bddeebe 100644
--- a/src/secondary_uncorrelated.cpp
+++ b/src/secondary_uncorrelated.cpp
@@ -52,7 +52,8 @@ UncorrelatedAngleEnergy::UncorrelatedAngleEnergy(hid_t group)
 }
 
 void
-UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
+UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
+  uint64_t * prn_seeds, int stream) const
 {
   // Sample cosine of scattering angle
   if (fission_) {
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 592878426..86f20a93d 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -476,7 +476,7 @@ void initialize_history(Particle* p, int64_t index_source)
   // set random number seed
   int64_t particle_seed = (simulation::total_gen + overall_generation() - 1)
     * settings::n_particles + p->id_;
-  set_particle_seed(particle_seed, p->prn_seeds);
+  set_particle_seed(particle_seed, p->prn_seeds_);
 
   // set particle trace
   simulation::trace = false;