diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h
index 9c04ff8c70..2d1ebc8f3b 100644
--- a/include/openmc/angle_energy.h
+++ b/include/openmc/angle_energy.h
@@ -1,6 +1,8 @@
 #ifndef OPENMC_ANGLE_ENERGY_H
 #define OPENMC_ANGLE_ENERGY_H
 
+#include <cstdint>
+
 namespace openmc {
 
 //==============================================================================
@@ -12,7 +14,8 @@ namespace openmc {
 
 class AngleEnergy {
 public:
-  virtual void sample(double E_in, double& E_out, double& mu) const = 0;
+  virtual void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+    int stream) const = 0;
   virtual ~AngleEnergy() = default;
 };
 
diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h
index 81ea039950..8ff5b0fc04 100644
--- a/include/openmc/distribution_energy.h
+++ b/include/openmc/distribution_energy.h
@@ -36,6 +36,8 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] prn_seeds Pseudorandom number genererator seed array
+  //! \param[in] stream Pseudorandom number genererator stream index
   //! \return Sampled energy in [eV]
   double sample(double E, uint64_t * prn_seeds, int stream) const;
 private:
@@ -55,6 +57,8 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] prn_seeds Pseudorandom number genererator seed array
+  //! \param[in] stream Pseudorandom number genererator stream index
   //! \return Sampled energy in [eV]
   double sample(double E, uint64_t * prn_seeds, int stream) const;
 private:
@@ -74,6 +78,8 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] prn_seeds Pseudorandom number genererator seed array
+  //! \param[in] stream Pseudorandom number genererator stream index
   //! \return Sampled energy in [eV]
   double sample(double E, uint64_t * prn_seeds, int stream) const;
 private:
@@ -103,6 +109,8 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] prn_seeds Pseudorandom number genererator seed array
+  //! \param[in] stream Pseudorandom number genererator stream index
   //! \return Sampled energy in [eV]
   double sample(double E, uint64_t * prn_seeds, int stream) const;
 private:
@@ -121,6 +129,8 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] prn_seeds Pseudorandom number genererator seed array
+  //! \param[in] stream Pseudorandom number genererator stream index
   //! \return Sampled energy in [eV]
   double sample(double E, uint64_t * prn_seeds, int stream) const;
 private:
@@ -139,6 +149,8 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] prn_seeds Pseudorandom number genererator seed array
+  //! \param[in] stream Pseudorandom number genererator stream index
   //! \return Sampled energy in [eV]
   double sample(double E, uint64_t * prn_seeds, int stream) const;
 private:
diff --git a/include/openmc/physics.h b/include/openmc/physics.h
index 8dc2650aee..35063ca0e6 100644
--- a/include/openmc/physics.h
+++ b/include/openmc/physics.h
@@ -44,7 +44,7 @@ void sample_positron_reaction(Particle* p);
 //!
 //! \param[in] p Particle
 //! \return Index in the data::nuclides vector
-int sample_nuclide(const Particle* p);
+int sample_nuclide(Particle* p);
 
 //! Determine the average total, prompt, and delayed neutrons produced from
 //! fission and creates appropriate bank sites.
@@ -53,9 +53,9 @@ void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 
 int sample_element(Particle* p);
 
-Reaction* sample_fission(int i_nuclide, const Particle* p);
+Reaction* sample_fission(int i_nuclide, Particle* p);
 
-void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_product);
+void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product);
 
 void absorption(Particle* p, int i_nuclide);
 
diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h
index a0daf803c2..4a6e13a836 100644
--- a/include/openmc/random_lcg.h
+++ b/include/openmc/random_lcg.h
@@ -37,7 +37,7 @@ extern "C" double prn(uint64_t * seeds, int stream);
 //==============================================================================
 
 //extern "C" double future_prn(int64_t n);
-extern "C" double future_prn(int64_t n, uint64_t seed);
+extern "C" double future_prn(int64_t n, uint64_t * prn_seeds, int stream);
 
 //==============================================================================
 //! Set the RNG seed to a unique value based on the ID of the particle.
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index 8d93c21f52..b17cbca564 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -103,7 +103,7 @@ public:
   //! \param[out] elastic Thermal elastic scattering cross section
   //! \param[out] inelastic Thermal inelastic scattering cross section
   void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic,
-                    double* inelastic) const;
+                    double* inelastic, uint64_t * prn_seeds, int stream) const;
 
   //! Determine whether table applies to a particular nuclide
   //!
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 1db293b0a0..48b86b8cc7 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -25,7 +25,7 @@ DiscretePhoton::DiscretePhoton(hid_t group)
   read_attribute(group, "atomic_weight_ratio", A_);
 }
 
-double DiscretePhoton::sample(double E) const
+double DiscretePhoton::sample(double E, uint64_t * prn_seeds, int stream) const
 {
   if (primary_flag_ == 2) {
     return energy_ + A_/(A_+ 1)*E;
@@ -44,7 +44,7 @@ LevelInelastic::LevelInelastic(hid_t group)
   read_attribute(group, "mass_ratio", mass_ratio_);
 }
 
-double LevelInelastic::sample(double E) const
+double LevelInelastic::sample(double E, uint64_t * prn_seeds, int stream ) const
 {
   return mass_ratio_*(E - threshold_);
 }
diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp
index 4468febd33..df2f8d682d 100644
--- a/src/distribution_multi.cpp
+++ b/src/distribution_multi.cpp
@@ -67,7 +67,7 @@ Direction PolarAzimuthal::sample(uint64_t * prn_seeds, int stream) const
   // TODO: apply this change directly to rotate_angle
   if (u_ref_.x == 0 && u_ref_.y == 0) phi += 0.5*PI;
 
-  return rotate_angle(u_ref_, mu, &phi);
+  return rotate_angle(u_ref_, mu, &phi, prn_seeds, stream);
 }
 
 //==============================================================================
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 1066b900c7..ccfd211896 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -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);
+  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);
diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp
index ce65556084..27ba1784dd 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);
+  set_particle_seed(particle_seed, p.prn_seeds);
 
   // Transport neutron
   p.transport();
diff --git a/src/physics.cpp b/src/physics.cpp
index 165975ea78..a17c0a8a43 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -398,7 +398,7 @@ void sample_positron_reaction(Particle* p)
   p->event_ = EVENT_ABSORB;
 }
 
-int sample_nuclide(const Particle* p)
+int sample_nuclide(Particle* p)
 {
   // Sample cumulative distribution function
   double cutoff = prn(p->prn_seeds, p->stream) * p->macro_xs_.total;
@@ -455,7 +455,7 @@ int sample_element(Particle* p)
   fatal_error("Did not sample any element during collision.");
 }
 
-Reaction* sample_fission(int i_nuclide, const Particle* p)
+Reaction* sample_fission(int i_nuclide, Particle* p)
 {
   // Get pointer to nuclide
   const auto& nuc {data::nuclides[i_nuclide]};
@@ -500,7 +500,7 @@ Reaction* sample_fission(int i_nuclide, const Particle* p)
   throw std::runtime_error{"No fission reaction was sampled for " + nuc->name_};
 }
 
-void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_product)
+void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product)
 {
   // Get grid index and interpolation factor and sample photon production cdf
   int i_temp = p->neutron_xs_[i_nuclide].index_temp;
diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp
index 62e8e17a35..c3b0c6a62f 100644
--- a/src/random_lcg.cpp
+++ b/src/random_lcg.cpp
@@ -7,7 +7,6 @@ namespace openmc {
 
 
 // Constants
-extern "C" const int N_STREAMS         {6};
 extern "C" const int STREAM_TRACKING   {0};
 extern "C" const int STREAM_TALLIES    {1};
 extern "C" const int STREAM_SOURCE     {2};
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index 9f4f652a06..5d32d3e860 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -33,13 +33,13 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu,
   double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_);
 
   // x is essentially a Maxwellian distribution
-  double x = maxwell_spectrum(1.0);
+  double x = maxwell_spectrum(1.0, prn_seeds, stream);
 
   double y;
   double r1, r2, r3, r4, r5, r6;
   switch (n_bodies_) {
   case 3:
-    y = maxwell_spectrum(1.0);
+    y = maxwell_spectrum(1.0, prn_seeds, stream);
     break;
   case 4:
     r1 = prn(prn_seeds, stream);