diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
index 603d8c414..b240034d8 100644
--- a/include/openmc/secondary_correlated.h
+++ b/include/openmc/secondary_correlated.h
@@ -38,7 +38,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+		  int stream) const override;
 
   // energy property
   std::vector<double>& energy() { return energy_; }
diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h
index 0217c7e77..90ad1e710 100644
--- a/include/openmc/secondary_kalbach.h
+++ b/include/openmc/secondary_kalbach.h
@@ -29,7 +29,10 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] prn_seeds Pseudorandom stream seed array
+  //! \param[in] stream Pseudorandom stream index
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+		  int sample) const override;
 private:
   //! Outgoing energy/angle at a single incoming energy
   struct KMTable {
diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h
index eb90bd8a0..13a13cc19 100644
--- a/include/openmc/secondary_nbody.h
+++ b/include/openmc/secondary_nbody.h
@@ -24,7 +24,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+  int stream) const override;
 private:
   int n_bodies_; //!< Number of particles distributed
   double mass_ratio_; //!< Total mass of particles [neutron mass]
diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h
index 7b876334f..17b5fe9f3 100644
--- a/include/openmc/secondary_thermal.h
+++ b/include/openmc/secondary_thermal.h
@@ -31,7 +31,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+    int stream) const override;
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
 };
@@ -51,7 +52,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+    int stream) const override;
 private:
   double debye_waller_;
 };
@@ -72,7 +74,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+    int stream) const override;
 private:
   const std::vector<double>& energy_; //!< Energies at which cosines are tabulated
   xt::xtensor<double, 2> mu_out_; //!< Cosines for each incident energy
@@ -94,7 +97,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+    int stream) const override;
 private:
   const std::vector<double>& energy_; //!< Incident energies
   xt::xtensor<double, 2> energy_out_; //!< Outgoing energies for each incident energy
@@ -117,7 +121,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+    int stream) const override;
 private:
   //! Secondary energy/angle distribution
   struct DistEnergySab {
diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h
index e430c75c4..69d097af6 100644
--- a/include/openmc/secondary_uncorrelated.h
+++ b/include/openmc/secondary_uncorrelated.h
@@ -29,7 +29,8 @@ public:
   //! \param[in] E_in Incoming energy in [eV]
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[out] mu Outgoing cosine with respect to current direction
-  void sample(double E_in, double& E_out, double& mu) const override;
+  void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
+    int stream) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }
diff --git a/include/openmc/source.h b/include/openmc/source.h
index bdd72e5bf..639f6ecab 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -39,7 +39,7 @@ public:
 
   //! Sample from the external source distribution
   //! \return Sampled site
-  Particle::Bank sample() const;
+  Particle::Bank sample(uint64_t * prn_seed, int stream) const;
 
   // Properties
   double strength() const { return strength_; }
@@ -61,7 +61,7 @@ extern "C" void initialize_source();
 //! Sample a site from all external source distributions in proportion to their
 //! source strength
 //! \return Sampled source site
-Particle::Bank sample_external_source();
+Particle::Bank sample_external_source(uint64_t * prn_seeds);
 
 //! Fill source bank at end of generation for fixed source simulations
 void fill_source_bank_fixedsource();
diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 5272dab3f..8060afcd1 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -116,7 +116,8 @@ public:
   //! \return Outgoing direction of the ray
   virtual Direction reflect(Position r, Direction u) const;
 
-  virtual Direction diffuse_reflect(Position r, Direction u) const;
+  virtual Direction diffuse_reflect(Position r, Direction u, uint64_t * prn_seeds,
+    int stream) const;
 
   //! Evaluate the equation describing the surface.
   //!
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index cb6a7ce7b..8d93c21f5 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -64,8 +64,10 @@ public:
   //! \param[in] E_in Incident neutron energy in [eV]
   //! \param[out] E_out Outgoing neutron energy in [eV]
   //! \param[out] mu Outgoing scattering angle cosine
+  //! \param[inout] prn_seeds Pseudorandom seed array
+  //! \param[in] stream Pseudorandom stream index
   void sample(const NuclideMicroXS& micro_xs, double E_in,
-              double* E_out, double* mu);
+              double* E_out, double* mu, uint64_t * prn_seeds, int stream);
 private:
   struct Reaction {
     // Default constructor
diff --git a/src/particle.cpp b/src/particle.cpp
index ab058a6c3..10c58239f 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -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());
+      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/secondary_correlated.cpp b/src/secondary_correlated.cpp
index 0a5d8804d..4c3b93262 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -152,14 +152,15 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
   } // incoming energies
 }
 
-void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
+void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
+  uint64_t * prn_seeds, int stream) const
 {
   // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   // Before the secondary distribution refactor, an isotropic polar cosine was
   // always sampled but then overwritten with the polar cosine sampled from the
   // correlated distribution. To preserve the random number stream, we keep
   // this dummy sampling here but can remove it later (will change answers)
-  mu = 2.0*prn() - 1.0;
+  mu = 2.0*prn(prn_seeds, stream) - 1.0;
   // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
   // Find energy bin and calculate interpolation factor -- if the energy is
@@ -179,7 +180,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
   }
 
   // Sample between the ith and [i+1]th bin
-  int l = r > prn() ? i + 1 : i;
+  int l = r > prn(prn_seeds, stream) ? i + 1 : i;
 
   // Interpolation for energy E1 and EK
   int n_energy_out = distribution_[i].e_out.size();
@@ -198,7 +199,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
   n_discrete = distribution_[l].n_discrete;
-  double r1 = prn();
+  double r1 = prn(prn_seeds, stream);
   double c_k = distribution_[l].c[0];
   int k = 0;
   int end = n_energy_out - 2;
@@ -259,9 +260,9 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
 
   // Find correlated angular distribution for closest outgoing energy bin
   if (r1 - c_k < c_k1 - r1) {
-    mu = distribution_[l].angle[k]->sample();
+    mu = distribution_[l].angle[k]->sample(prn_seeds, stream);
   } else {
-    mu = distribution_[l].angle[k + 1]->sample();
+    mu = distribution_[l].angle[k + 1]->sample(prn_seeds, stream);
   }
 }
 
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index e3974f357..781923720 100644
--- a/src/secondary_kalbach.cpp
+++ b/src/secondary_kalbach.cpp
@@ -113,14 +113,14 @@ KalbachMann::KalbachMann(hid_t group)
   } // incoming energies
 }
 
-void KalbachMann::sample(double E_in, double& E_out, double& mu) const
+void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
 {
   // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   // Before the secondary distribution refactor, an isotropic polar cosine was
   // always sampled but then overwritten with the polar cosine sampled from the
   // correlated distribution. To preserve the random number stream, we keep
   // this dummy sampling here but can remove it later (will change answers)
-  mu = 2.0*prn() - 1.0;
+  mu = 2.0*prn(prn_seeds, stream) - 1.0;
   // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
   // Find energy bin and calculate interpolation factor -- if the energy is
@@ -140,7 +140,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
   }
 
   // Sample between the ith and [i+1]th bin
-  int l = r > prn() ? i + 1 : i;
+  int l = r > prn(prn_seeds, stream) ? i + 1 : i;
 
   // Interpolation for energy E1 and EK
   int n_energy_out = distribution_[i].e_out.size();
@@ -159,7 +159,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
   n_discrete = distribution_[l].n_discrete;
-  double r1 = prn();
+  double r1 = prn(prn_seeds, stream);
   double c_k = distribution_[l].c[0];
   int k = 0;
   int end = n_energy_out - 2;
@@ -229,11 +229,11 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu) const
   }
 
   // Sampled correlated angle from Kalbach-Mann parameters
-  if (prn() > km_r) {
-    double T = (2.0*prn() - 1.0) * std::sinh(km_a);
+  if (prn(prn_seeds, stream) > km_r) {
+    double T = (2.0*prn(prn_seeds, stream) - 1.0) * std::sinh(km_a);
     mu = std::log(T + std::sqrt(T*T + 1.0))/km_a;
   } else {
-    double r1 = prn();
+    double r1 = prn(prn_seeds, stream);
     mu = std::log(r1*std::exp(km_a) + (1.0 - r1)*std::exp(-km_a))/km_a;
   }
 }
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index 564fca73d..9f4f652a0 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -21,11 +21,12 @@ NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)
   read_attribute(group, "q_value", Q_);
 }
 
-void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const
+void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu,
+  uint64_t * prn_seeds, int stream) const
 {
   // By definition, the distribution of the angle is isotropic for an N-body
   // phase space distribution
-  mu = 2.0*prn() - 1.0;
+  mu = 2.0*prn(prn_seeds, stream) - 1.0;
 
   // Determine E_max parameter
   double Ap = mass_ratio_;
@@ -41,18 +42,18 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const
     y = maxwell_spectrum(1.0);
     break;
   case 4:
-    r1 = prn();
-    r2 = prn();
-    r3 = prn();
+    r1 = prn(prn_seeds, stream);
+    r2 = prn(prn_seeds, stream);
+    r3 = prn(prn_seeds, stream);
     y = -std::log(r1*r2*r3);
     break;
   case 5:
-    r1 = prn();
-    r2 = prn();
-    r3 = prn();
-    r4 = prn();
-    r5 = prn();
-    r6 = prn();
+    r1 = prn(prn_seeds, stream);
+    r2 = prn(prn_seeds, stream);
+    r3 = prn(prn_seeds, stream);
+    r4 = prn(prn_seeds, stream);
+    r5 = prn(prn_seeds, stream);
+    r6 = prn(prn_seeds, stream);
     y = -std::log(r1*r2*r3*r4) - std::log(r5) * std::pow(std::cos(PI/2.0*r6), 2);
     break;
   default:
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
index 8c9cc6967..7c3af2993 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -32,7 +32,7 @@ CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs)
 { }
 
 void
-CoherentElasticAE::sample(double E_in, double& E_out, double& mu) 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;
@@ -42,7 +42,7 @@ CoherentElasticAE::sample(double E_in, double& E_out, double& mu) const
 
   // Sample a Bragg edge between 1 and i
   const auto& factors = xs_.factors();
-  double prob = prn() * factors[i+1];
+  double prob = prn(prn_seeds, stream) * factors[i+1];
   int k = 0;
   if (prob >= factors.front()) {
     k = lower_bound_index(factors.begin(), factors.begin() + (i+1), prob);
@@ -65,11 +65,11 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group)
 }
 
 void
-IncoherentElasticAE::sample(double E_in, double& E_out, double& mu) 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_;
-  mu = std::log(1.0 + prn()*(std::exp(2.0*c) - 1))/c - 1.0;
+  mu = std::log(1.0 + prn(prn_seeds, stream)*(std::exp(2.0*c) - 1))/c - 1.0;
 
   // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
   E_out = E_in;
@@ -87,7 +87,7 @@ IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group,
 }
 
 void
-IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu) 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;
@@ -98,7 +98,7 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu) cons
   // incoming energies.
 
   // Sample outgoing cosine bin
-  int k = prn() * mu_out_.shape()[1];
+  int k = prn(prn_seeds, stream) * mu_out_.shape()[1];
 
   // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
   double mu_ik  = mu_out_(i, k);
@@ -125,7 +125,7 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group,
 }
 
 void
-IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu) 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;
@@ -143,10 +143,10 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu) co
   int n = energy_out_.shape()[1];
   if (!skewed_) {
     // All bins equally likely
-    j = prn() * n;
+    j = prn(prn_seeds, stream) * n;
   } else {
     // Distribution skewed away from edge points
-    double r = prn() * (n - 3);
+    double r = prn(prn_seeds, stream) * (n - 3);
     if (r > 1.0) {
       // equally likely N-4 middle bins
       j = r + 1;
@@ -174,7 +174,7 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu) co
 
   // Sample outgoing cosine bin
   int m = mu_out_.shape()[2];
-  int k = prn() * m;
+  int k = prn(prn_seeds, stream) * m;
 
   // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
   double mu_ijk  = mu_out_(i, j, k);
@@ -228,7 +228,7 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
 }
 
 void
-IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) 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;
@@ -236,7 +236,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const
   get_energy_index(energy_, E_in, i, f);
 
   // Sample between ith and [i+1]th bin
-  int l = f > prn() ? i + 1 : i;
+  int l = f > prn(prn_seeds, stream) ? i + 1 : i;
 
   // Determine endpoints on grid i
   auto n = distribution_[i].e_out.size();
@@ -254,7 +254,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const
   // Determine outgoing energy bin
   // (First reset n_energy_out to the right value)
   n = distribution_[l].n_e_out;
-  double r1 = prn();
+  double r1 = prn(prn_seeds, stream);
   double c_j = distribution_[l].e_out_cdf[0];
   double c_j1;
   std::size_t j;
@@ -293,7 +293,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const
 
   // Sample outgoing cosine bin
   int n_mu = distribution_[l].mu.shape()[1];
-  std::size_t k = prn() * n_mu;
+  std::size_t k = prn(prn_seeds, stream) * n_mu;
 
   // Rather than use the sampled discrete mu directly, it is smeared over
   // a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
@@ -321,7 +321,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const
   }
 
   // Smear angle
-  mu += std::min(mu - mu_left, mu_right - mu)*(prn() - 0.5);
+  mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seeds, stream) - 0.5);
 }
 
 } // namespace openmc
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
index 3192c1e04..1db62bf2f 100644
--- a/src/secondary_uncorrelated.cpp
+++ b/src/secondary_uncorrelated.cpp
@@ -52,7 +52,7 @@ UncorrelatedAngleEnergy::UncorrelatedAngleEnergy(hid_t group)
 }
 
 void
-UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) 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_) {
@@ -61,14 +61,14 @@ UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
     mu = 1.0;
     // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   } else if (!angle_.empty()) {
-    mu = angle_.sample(E_in);
+    mu = angle_.sample(E_in, prn_seeds, stream);
   } else {
     // no angle distribution given => assume isotropic for all energies
-    mu = 2.0*prn() - 1.0;
+    mu = 2.0*prn(prn_seeds, stream) - 1.0;
   }
 
   // Sample outgoing energy
-  E_out = energy_->sample(E_in);
+  E_out = energy_->sample(E_in, prn_seeds, stream);
 }
 
 } // namespace openmc
diff --git a/src/source.cpp b/src/source.cpp
index 251e295ec..860584a2c 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -142,7 +142,7 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
 }
 
 
-Particle::Bank SourceDistribution::sample() const
+Particle::Bank SourceDistribution::sample(uint64_t * prn_seeds, int stream) const
 {
   Particle::Bank site;
 
@@ -158,7 +158,7 @@ Particle::Bank SourceDistribution::sample() const
     site.particle = particle_;
 
     // Sample spatial distribution
-    site.r = space_->sample();
+    site.r = space_->sample(prn_seeds, stream);
     double xyz[] {site.r.x, site.r.y, site.r.z};
 
     // Now search to see if location exists in geometry
@@ -200,7 +200,7 @@ Particle::Bank SourceDistribution::sample() const
   ++n_accept;
 
   // Sample angle
-  site.u = angle_->sample();
+  site.u = angle_->sample(prn_seeds, stream);
 
   // Check for monoenergetic source above maximum particle energy
   auto p = static_cast<int>(particle_);
@@ -218,7 +218,7 @@ Particle::Bank SourceDistribution::sample() const
 
   while (true) {
     // Sample energy spectrum
-    site.E = energy_->sample();
+    site.E = energy_->sample(prn_seeds, stream);
 
     // Resample if energy falls outside minimum or maximum particle energy
     if (site.E < data::energy_max[p] && site.E > data::energy_min[p]) break;
@@ -270,10 +270,11 @@ void initialize_source()
       // initialize random number seed
       int64_t id = simulation::total_gen*settings::n_particles +
         simulation::work_index[mpi::rank] + i + 1;
-      set_particle_seed(id);
+	  uint64_t prn_seeds[N_STREAMS];
+      set_particle_seed(id, prn_seeds);
 
       // sample external source distribution
-      simulation::source_bank[i] = sample_external_source();
+      simulation::source_bank[i] = sample_external_source(prn_seeds);
     }
   }
 
@@ -287,10 +288,10 @@ void initialize_source()
   }
 }
 
-Particle::Bank sample_external_source()
+Particle::Bank sample_external_source(uint64_t * prn_seeds)
 {
   // Set the random number generator to the source stream.
-  prn_set_stream(STREAM_SOURCE);
+  int stream = STREAM_SOURCE;
 
   // Determine total source strength
   double total_strength = 0.0;
@@ -300,7 +301,7 @@ Particle::Bank sample_external_source()
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    double xi = prn()*total_strength;
+    double xi = prn(prn_seeds, stream)*total_strength;
     double c = 0.0;
     for (; i < model::external_sources.size(); ++i) {
       c += model::external_sources[i].strength();
@@ -309,7 +310,7 @@ Particle::Bank sample_external_source()
   }
 
   // Sample source site from i-th source distribution
-  Particle::Bank site {model::external_sources[i].sample()};
+  Particle::Bank site {model::external_sources[i].sample(prn_seeds, stream)};
 
   // If running in MG, convert site.E to group
   if (!settings::run_CE) {
@@ -318,9 +319,6 @@ Particle::Bank sample_external_source()
     site.E = data::mg.num_energy_groups_ - site.E - 1.;
   }
 
-  // Set the random number generator back to the tracking stream.
-  prn_set_stream(STREAM_TRACKING);
-
   return site;
 }
 
@@ -336,10 +334,11 @@ 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;
-      set_particle_seed(id);
+	  uint64_t prn_seeds[N_STREAMS];
+      set_particle_seed(id, prn_seeds);
 
       // sample external source distribution
-      simulation::source_bank[i] = sample_external_source();
+      simulation::source_bank[i] = sample_external_source(prn_seeds);
     }
   }
 }
diff --git a/src/surface.cpp b/src/surface.cpp
index 235465952..7c18d1c34 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -197,7 +197,7 @@ Surface::reflect(Position r, Direction u) const
 }
 
 Direction
-Surface::diffuse_reflect(Position r, Direction u) const
+Surface::diffuse_reflect(Position r, Direction u, uint64_t * prn_seeds, int stream) const
 {
   // Diffuse reflect direction according to the normal.
   // cosine distribution 
@@ -208,10 +208,10 @@ Surface::diffuse_reflect(Position r, Direction u) const
   
   // sample from inverse function, u=sqrt(rand) since p(u)=2u, so F(u)=u^2 
   const double mu = (projection>=0.0) ? 
-                  -std::sqrt(prn()) : std::sqrt(prn());  
+                  -std::sqrt(prn(prn_seeds, stream)) : std::sqrt(prn(prn_seeds, stream));  
   
   // sample azimuthal distribution uniformly 
-  u = rotate_angle(n, mu, nullptr);
+  u = rotate_angle(n, mu, nullptr, prn_seeds, stream);
   
   // normalize the direction 
   return u/u.norm();    
diff --git a/src/thermal.cpp b/src/thermal.cpp
index 49b522376..0af5b1614 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -150,7 +150,8 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& tem
 
 void
 ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
-                                double* elastic, double* inelastic) const
+                                double* elastic, double* inelastic,
+                                uint64_t * prn_seeds, int stream) const
 {
   // Determine temperature for S(a,b) table
   double kT = sqrtkT*sqrtkT;
@@ -172,7 +173,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
 
     // Randomly sample between temperature i and i+1
     double f = (kT - kTs_[i]) / (kTs_[i+1] - kTs_[i]);
-    if (f > prn()) ++i;
+    if (f > prn(prn_seeds, stream)) ++i;
   }
 
   // Set temperature index
@@ -265,13 +266,13 @@ ThermalData::calculate_xs(double E, double* elastic, double* inelastic) const
 
 void
 ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
-                    double* E_out, double* mu)
+                    double* E_out, double* mu, uint64_t * prn_seeds, int stream)
 {
   // Determine whether inelastic or elastic scattering will occur
-  if (prn() < micro_xs.thermal_elastic / micro_xs.thermal) {
-    elastic_.distribution->sample(E, *E_out, *mu);
+  if (prn(prn_seeds, stream) < micro_xs.thermal_elastic / micro_xs.thermal) {
+    elastic_.distribution->sample(E, *E_out, *mu, prn_seeds, stream);
   } else {
-    inelastic_.distribution->sample(E, *E_out, *mu);
+    inelastic_.distribution->sample(E, *E_out, *mu, prn_seeds, stream);
   }
 
   // Because of floating-point roundoff, it may be possible for mu to be
diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp
index 419bf0ea5..0eb2df6f4 100644
--- a/src/volume_calc.cpp
+++ b/src/volume_calc.cpp
@@ -125,15 +125,16 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
       std::vector<std::vector<int>> hits(n);
       Particle p;
 
-      prn_set_stream(STREAM_VOLUME);
+      uint64_t prn_seeds[N_STREAMS];
+      int stream = STREAM_VOLUME;
 
       // Sample locations and count hits
       #pragma omp for
       for (size_t i = i_start; i < i_end; i++) {
-        set_particle_seed(iterations * n_samples_ + i);
+        set_particle_seed(iterations * n_samples_ + i, prn_seeds);
 
         p.n_coord_ = 1;
-        Position xi {prn(), prn(), prn()};
+        Position xi {prn(prn_seeds, stream), prn(prn_seeds, stream), prn(prn_seeds, stream)};
         p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
         p.u() = {0.5, 0.5, 0.5};
 
@@ -203,7 +204,6 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
           }
         }
       }
-      prn_set_stream(STREAM_TRACKING);
     } // omp parallel
 
     // Reduce hits onto master process