diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h
index c674035bf..c88acd569 100644
--- a/include/openmc/angle_energy.h
+++ b/include/openmc/angle_energy.h
@@ -14,8 +14,8 @@ namespace openmc {
 
 class AngleEnergy {
 public:
-  virtual void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const = 0;
+  virtual void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) const = 0;
   virtual ~AngleEnergy() = default;
 };
 
diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 8d5516ddf..d62ba1624 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -21,7 +21,7 @@ namespace openmc {
 class Distribution {
 public:
   virtual ~Distribution() = default;
-  virtual double sample(uint64_t* prn_seeds, int stream) const = 0;
+  virtual double sample(uint64_t* prn_seed) const = 0;
 };
 
 //==============================================================================
@@ -34,10 +34,9 @@ public:
   Discrete(const double* x, const double* p, int n);
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 
   // Properties
   const std::vector<double>& x() const { return x_; }
@@ -60,10 +59,9 @@ public:
   Uniform(double a, double b) : a_{a}, b_{b} {};
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 private:
   double a_; //!< Lower bound of distribution
   double b_; //!< Upper bound of distribution
@@ -79,10 +77,9 @@ public:
   Maxwell(double theta) : theta_{theta} { };
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 private:
   double theta_; //!< Factor in exponential [eV]
 };
@@ -97,10 +94,9 @@ public:
   Watt(double a, double b) : a_{a}, b_{b} { };
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 private:
   double a_; //!< Factor in exponential [eV]
   double b_; //!< Factor in square root [1/eV]
@@ -116,10 +112,9 @@ public:
   Normal(double mean_value, double std_dev) : mean_value_{mean_value}, std_dev_{std_dev} { };
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 private:
   double mean_value_;    //!< middle of distribution [eV]
   double std_dev_; //!< standard deviation [eV]
@@ -136,10 +131,9 @@ public:
   Muir(double e0, double m_rat, double kt) : e0_{e0}, m_rat_{m_rat}, kt_{kt} { };
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 private:
   // example DT fusion m_rat = 5 (D = 2 + T = 3)
   // ion temp = 20000 eV
@@ -160,10 +154,9 @@ public:
           const double* c=nullptr);
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 
   // x property
   std::vector<double>& x() { return x_; }
@@ -192,10 +185,9 @@ public:
   Equiprobable(const double* x, int n) : x_{x, x+n} { };
 
   //! Sample a value from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample(uint64_t* prn_seeds, int stream) const;
+  double sample(uint64_t* prn_seed) const;
 private:
   std::vector<double> x_; //! Possible outcomes
 };
diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h
index 307dca8d6..23d7aa008 100644
--- a/include/openmc/distribution_angle.h
+++ b/include/openmc/distribution_angle.h
@@ -23,10 +23,9 @@ public:
 
   //! Sample an angle given an incident particle energy
   //! \param[in] E Particle energy in [eV]
-  //! \param[inout] prn_seeds Array of pseudorandom number seeds
-  //! \param[in] stream Pseudorandom stream index
+  //! \param[inout] prn_seed pseudorandom number seed pointer
   //! \return Cosine of the angle in the range [-1,1]
-  double sample(double E, uint64_t* prn_seeds, int stream) const;
+  double sample(double E, uint64_t* prn_seed) const;
 
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty
diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h
index fe4524690..dd59a2359 100644
--- a/include/openmc/distribution_energy.h
+++ b/include/openmc/distribution_energy.h
@@ -22,7 +22,7 @@ namespace openmc {
 
 class EnergyDistribution {
 public:
-  virtual double sample(double E, uint64_t* prn_seeds, int stream) const = 0;
+  virtual double sample(double E, uint64_t* prn_seed) const = 0;
   virtual ~EnergyDistribution() = default;
 };
 
@@ -36,10 +36,9 @@ 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
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E, uint64_t* prn_seeds, int stream) const;
+  double sample(double E, uint64_t* prn_seed) const;
 private:
   int primary_flag_; //!< Indicator of whether the photon is a primary or
                      //!< non-primary photon.
@@ -57,10 +56,9 @@ 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
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E, uint64_t* prn_seeds, int stream) const;
+  double sample(double E, uint64_t* prn_seed) const;
 private:
   double threshold_; //!< Energy threshold in lab, (A + 1)/A * |Q|
   double mass_ratio_; //!< (A/(A+1))^2
@@ -78,10 +76,9 @@ 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
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E, uint64_t* prn_seeds, int stream) const;
+  double sample(double E, uint64_t* prn_seed) const;
 private:
   //! Outgoing energy for a single incoming energy
   struct CTTable {
@@ -109,10 +106,9 @@ 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
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E, uint64_t* prn_seeds, int stream) const;
+  double sample(double E, uint64_t* prn_seed) const;
 private:
   Tabulated1D theta_; //!< Incoming energy dependent parameter
   double u_; //!< Restriction energy
@@ -129,10 +125,9 @@ 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
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E, uint64_t* prn_seeds, int stream) const;
+  double sample(double E, uint64_t* prn_seed) const;
 private:
   Tabulated1D theta_; //!< Incoming energy dependent parameter
   double u_; //!< Restriction energy
@@ -149,10 +144,9 @@ 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
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E, uint64_t* prn_seeds, int stream) const;
+  double sample(double E, uint64_t* prn_seed) const;
 private:
   Tabulated1D a_; //!< Energy-dependent 'a' parameter
   Tabulated1D b_; //!< Energy-dependent 'b' parameter
diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h
index b8dcce267..08f1cd0b5 100644
--- a/include/openmc/distribution_multi.h
+++ b/include/openmc/distribution_multi.h
@@ -23,10 +23,9 @@ public:
   virtual ~UnitSphereDistribution() = default;
 
   //! Sample a direction from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Direction sampled
-  virtual Direction sample(uint64_t* prn_seeds, int stream) const = 0;
+  virtual Direction sample(uint64_t* prn_seed) const = 0;
 
   Direction u_ref_ {0.0, 0.0, 1.0};  //!< reference direction
 };
@@ -41,10 +40,9 @@ public:
   explicit PolarAzimuthal(pugi::xml_node node);
 
   //! Sample a direction from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Direction sampled
-  Direction sample(uint64_t* prn_seeds, int stream) const;
+  Direction sample(uint64_t* prn_seed) const;
 private:
   UPtrDist mu_;  //!< Distribution of polar angle
   UPtrDist phi_; //!< Distribution of azimuthal angle
@@ -59,10 +57,9 @@ public:
   Isotropic() { };
 
   //! Sample a direction from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled direction
-  Direction sample(uint64_t* prn_seeds, int stream) const;
+  Direction sample(uint64_t* prn_seed) const;
 };
 
 //==============================================================================
@@ -75,10 +72,9 @@ public:
   explicit Monodirectional(pugi::xml_node node) : UnitSphereDistribution{node} { };
 
   //! Sample a direction from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled direction
-  Direction sample(uint64_t* prn_seeds, int stream) const;
+  Direction sample(uint64_t* prn_seed) const;
 };
 
 using UPtrAngle = std::unique_ptr<UnitSphereDistribution>;
diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h
index c6cd1a26e..4895dd6ae 100644
--- a/include/openmc/distribution_spatial.h
+++ b/include/openmc/distribution_spatial.h
@@ -17,7 +17,7 @@ public:
   virtual ~SpatialDistribution() = default;
 
   //! Sample a position from the distribution
-  virtual Position sample(uint64_t* prn_seeds, int stream) const = 0;
+  virtual Position sample(uint64_t* prn_seed) const = 0;
 };
 
 //==============================================================================
@@ -29,10 +29,9 @@ public:
   explicit CartesianIndependent(pugi::xml_node node);
 
   //! Sample a position from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample(uint64_t* prn_seeds, int stream) const;
+  Position sample(uint64_t* prn_seed) const;
 private:
   UPtrDist x_; //!< Distribution of x coordinates
   UPtrDist y_; //!< Distribution of y coordinates
@@ -48,10 +47,9 @@ public:
   explicit SphericalIndependent(pugi::xml_node node);
 
   //! Sample a position from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample(uint64_t* prn_seeds, int stream) const;
+  Position sample(uint64_t* prn_seed) const;
 private:
   UPtrDist r_; //!< Distribution of r coordinates
   UPtrDist theta_; //!< Distribution of theta coordinates
@@ -68,10 +66,9 @@ public:
   explicit SpatialBox(pugi::xml_node node, bool fission=false);
 
   //! Sample a position from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample(uint64_t* prn_seeds, int stream) const;
+  Position sample(uint64_t* prn_seed) const;
 
   // Properties
   bool only_fissionable() const { return only_fissionable_; }
@@ -92,10 +89,9 @@ public:
   explicit SpatialPoint(pugi::xml_node node);
 
   //! Sample a position from the distribution
-  //! \param prn_seeds Array of pseudorandom number seeds
-  //! \param stream Pseudorandom stream index
+  //! \param prn_seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample(uint64_t* prn_seeds, int stream) const;
+  Position sample(uint64_t* prn_seed) const;
 private:
   Position r_; //!< Single position at which sites are generated
 };
diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index 5c7d58d58..1440d13b9 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -129,15 +129,14 @@ extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]);
 //! \param mu     The cosine of angle in lab or CM
 //! \param phi    The azimuthal angle; will randomly chosen angle if a nullptr
 //!   is passed
-//! \param prn_seeds A pointer to the array of pseudorandom seeds
-//! \param stream The pseudorandom stream index with which to sample from
+//! \param prn_seed A pointer to the pseudorandom seed
 //==============================================================================
 
 extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi,
-  uint64_t* prn_seeds, int stream);
+  uint64_t* prn_seed);
 
 Direction rotate_angle(Direction u, double mu, const double* phi,
-  uint64_t* prn_streams, int stream);
+  uint64_t* prn_seed);
 
 //==============================================================================
 //! Samples an energy from the Maxwell fission distribution based on a direct
@@ -148,12 +147,11 @@ Direction rotate_angle(Direction u, double mu, const double* phi,
 //! rule C64 in the Monte Carlo Sampler LA-9721-MS.
 //!
 //! \param T The tabulated function of the incoming energy
-//! \param prn_seeds A pointer to the array of pseudorandom seeds
-//! \param stream The pseudorandom stream index with which to sample from
+//! \param prn_seed A pointer to the pseudorandom seed
 //! \return The sampled outgoing energy
 //==============================================================================
 
-extern "C" double maxwell_spectrum(double T, uint64_t* prn_seeds, int stream);
+extern "C" double maxwell_spectrum(double T, uint64_t* prn_seed);
 
 //==============================================================================
 //! Samples an energy from a Watt energy-dependent fission distribution.
@@ -165,12 +163,11 @@ extern "C" double maxwell_spectrum(double T, uint64_t* prn_seeds, int stream);
 //!
 //! \param a Watt parameter a
 //! \param b Watt parameter b
-//! \param prn_seeds A pointer to the array of pseudorandom seeds
-//! \param stream The pseudorandom stream index with which to sample from
+//! \param prn_seed A pointer to the pseudorandom seed
 //! \return The sampled outgoing energy
 //==============================================================================
 
-extern "C" double watt_spectrum(double a, double b, uint64_t* prn_seeds, int stream);
+extern "C" double watt_spectrum(double a, double b, uint64_t* prn_seed);
 
 //==============================================================================
 //! Samples an energy from the Gaussian energy-dependent fission distribution.
@@ -183,13 +180,11 @@ extern "C" double watt_spectrum(double a, double b, uint64_t* prn_seeds, int str
 //!
 //! @param mean mean of the Gaussian distribution
 //! @param std_dev standard deviation of the Gaussian distribution
-//! @param prn_seeds A pointer to the array of pseudorandom seeds
-//! @param stream The pseudorandom stream index with which to sample from
+//! @param prn_seed A pointer to the pseudorandom seed
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double normal_variate(double mean, double std_dev, uint64_t* prn_seeds,
-  int stream);
+extern "C" double normal_variate(double mean, double std_dev, uint64_t* prn_seed);
 
 //==============================================================================
 //! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
@@ -201,13 +196,12 @@ extern "C" double normal_variate(double mean, double std_dev, uint64_t* prn_seed
 //! @param e0 peak neutron energy [eV]
 //! @param m_rat ratio of the fusion reactants to AMU
 //! @param kt the ion temperature of the reactants [eV]
-//! @param prn_seeds A pointer to the array of pseudorandom seeds
-//! @param stream The pseudorandom stream index with which to sample from
+//! @param prn_seed A pointer to the pseudorandom seed
 //! @result The sampled outgoing energy
 //==============================================================================
 
 extern "C" double muir_spectrum(double e0, double m_rat, double kt,
-  uint64_t* prn_seeds, int stream);
+  uint64_t* prn_seed);
 
 //==============================================================================
 //! Doppler broadens the windowed multipole curvefit.
diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h
index f32415a24..131d248e1 100644
--- a/include/openmc/mgxs.h
+++ b/include/openmc/mgxs.h
@@ -154,10 +154,9 @@ class Mgxs {
     //! @param gin Incoming energy group.
     //! @param dg Sampled delayed group index.
     //! @param gout Sampled outgoing energy group.
-    //! @param prn_seeds Pseudorandom seed array.
-    //! @param stream Pseudorandom stream index.
+    //! @param prn_seed Pseudorandom seed pointer
     void
-    sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seeds, int stream);
+    sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed);
 
     //! \brief Samples the outgoing energy and angle from a scatter event.
     //!
@@ -165,11 +164,9 @@ class Mgxs {
     //! @param gout Sampled outgoing energy group.
     //! @param mu Sampled cosine of the change-in-angle.
     //! @param wgt Weight of the particle to be adjusted.
-    //! @param prn_seeds Array of pseudorandom stream seeds
-    //! @param stream Pseudorandom stream index
+    //! @param prn_seed Pseudorandom seed pointer.
     void
-    sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seeds,
-      int stream);
+    sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
 
     //! \brief Calculates cross section quantities needed for tracking.
     //!
diff --git a/include/openmc/photon.h b/include/openmc/photon.h
index 81def193b..c7bbeddd2 100644
--- a/include/openmc/photon.h
+++ b/include/openmc/photon.h
@@ -45,12 +45,12 @@ public:
   void calculate_xs(Particle& p) const;
 
   void compton_scatter(double alpha, bool doppler, double* alpha_out,
-    double* mu, int* i_shell, uint64_t* prn_seeds, int stream) const;
+    double* mu, int* i_shell, uint64_t* prn_seed) const;
 
-  double rayleigh_scatter(double alpha, uint64_t* prn_seeds, int stream) const;
+  double rayleigh_scatter(double alpha, uint64_t* prn_seed) const;
 
   void pair_production(double alpha, double* E_electron, double* E_positron,
-    double* mu_electron, double* mu_positron, uint64_t* prn_seeds, int stream) const;
+    double* mu_electron, double* mu_positron, uint64_t* prn_seed) const;
 
   void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const;
 
@@ -97,14 +97,14 @@ public:
 
 private:
   void compton_doppler(double alpha, double mu, double* E_out, int* i_shell,
-                       uint64_t* prn_seeds, int stream) const;
+                       uint64_t* prn_seed) const;
 };
 
 //==============================================================================
 // Non-member functions
 //==============================================================================
 
-std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seeds, int stream);
+std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seed);
 
 void free_memory_photon();
 
diff --git a/include/openmc/physics.h b/include/openmc/physics.h
index de77b046d..e4a6183be 100644
--- a/include/openmc/physics.h
+++ b/include/openmc/physics.h
@@ -72,17 +72,17 @@ void sab_scatter(int i_nuclide, int i_sab, Particle* p);
 //! dependence of cross sections in treating resonance elastic scattering such
 //! as the DBRC and a new, accelerated scheme are also implemented here.
 Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
-   Direction v_neut, double xs_eff, double kT, uint64_t* prn_seeds, int stream);
+   Direction v_neut, double xs_eff, double kT, uint64_t* prn_seed);
 
 //! samples a target velocity based on the free gas scattering formulation, used
 //! by most Monte Carlo codes, in which cross section is assumed to be constant
 //! in energy. Excellent documentation for this method can be found in
 //! FRA-TM-123.
 Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT,
-		uint64_t* prn_seeds, int stream);
+		uint64_t* prn_seed);
 
 void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in,
-  Particle::Bank* site, uint64_t* prn_seeds, int stream);
+  Particle::Bank* site, uint64_t* prn_seed);
 
 //! handles all reactions with a single secondary neutron (other than fission),
 //! i.e. level scattering, (n,np), (n,na), etc.
diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h
index 691bd24c5..53cb39f1e 100644
--- a/include/openmc/random_lcg.h
+++ b/include/openmc/random_lcg.h
@@ -21,42 +21,53 @@ constexpr int64_t DEFAULT_SEED = 1;
 
 //==============================================================================
 //! Generate a pseudo-random number using a linear congruential generator.
-//! @param prn_seeds Pseudorandom number seed array
-//! @param stream Pseudorandom number stream index
+//! @param prn_seed Pseudorandom number seed pointer
 //! @return A random number between 0 and 1
 //==============================================================================
 
-extern "C" double prn(uint64_t* seeds, int stream);
+extern "C" double prn(uint64_t* prn_seed);
 
 //==============================================================================
 //! Generate a random number which is 'n' times ahead from the current seed.
 //!
 //! The result of this function will be the same as the result from calling
-//! `prn()` 'n' times.
+//! `prn()` 'n' times, though without the side effect of altering the RNG
+//! state.
 //! @param n The number of RNG seeds to skip ahead by
-//! @param prn_seeds Pseudorandom number seed array
-//! @param stream Pseudorandom number stream index
+//! @param prn_seed Pseudorandom number seed
 //! @return A random number between 0 and 1
 //==============================================================================
 
-extern "C" double future_prn(int64_t n, uint64_t* prn_seeds, int stream);
+extern "C" double future_prn(int64_t n, uint64_t prn_seed);
 
 //==============================================================================
-//! Set the RNG seeds to unique values based on the ID of the particle.
+//! 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
 //==============================================================================
 
-extern "C" void set_particle_seed(int64_t id, uint64_t* prn_seeds );
+extern "C" void init_seed(int64_t id, uint64_t* prn_seeds, int offset );
 
 //==============================================================================
-//! Advance the random number seed 'n' times from the current seed.
+//! Set the RNG seeds to unique values based on the ID of the particle. This
+//! function initializes the seeds for all RNG streams of the particle via
+//! striding.
 //! @param prn_seeds Pseudorandom number seed array
-//! @param stream Pseudorandom number stream index
+//! @param id The particle ID
+//==============================================================================
+
+extern "C" void init_particle_seeds(int64_t id, uint64_t* prn_seeds );
+
+//==============================================================================
+//! Advance the random number seed 'n' times from the current seed. This
+//! differs from the future_prn() function in that this function does alter
+//! the RNG state.
+//! @param prn_seed Pseudorandom number seed pointer
 //! @param n The number of RNG seeds to skip ahead by
 //==============================================================================
 
-extern "C" void advance_prn_seed(int64_t n, uint64_t* prn_seeds, int stream);
+extern "C" void advance_prn_seed(int64_t n, uint64_t* prn_seed);
 
 //==============================================================================
 //! Advance a random number seed 'n' times.
@@ -68,7 +79,7 @@ extern "C" void advance_prn_seed(int64_t n, uint64_t* prn_seeds, int stream);
 //! @param seed The starting to seed to advance from
 //==============================================================================
 
-uint64_t future_seed(uint64_t n, uint64_t seed);
+uint64_t future_seed(uint64_t n, uint64_t prn_seed);
 
 //==============================================================================
 //                               API FUNCTIONS
diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h
index 48f7f5e10..c61eec6b0 100644
--- a/include/openmc/reaction_product.h
+++ b/include/openmc/reaction_product.h
@@ -42,10 +42,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
-  //! \param[inout] prn_seeds Pseudorandom array of stream seeds
-  //! \param[in] stream Psuedorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const;
+  //! \param[inout] prn_seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seed) const;
 
   Particle::Type particle_; //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h
index fb5d4cd7a..8f3a79b77 100644
--- a/include/openmc/scattdata.h
+++ b/include/openmc/scattdata.h
@@ -65,11 +65,9 @@ class ScattData {
     //! @param gout Sampled outgoing energy group.
     //! @param mu Sampled cosine of the change-in-angle.
     //! @param wgt Weight of the particle to be adjusted.
-    //! @param prn_seeds Array of pseudorandom stream seeds
-    //! @param stream Pseudorandom stream index
+    //! @param prn_seed Pseudorandom number seed pointer
     virtual void
-    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seeds,
-      int stream) = 0;
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed) = 0;
 
     //! \brief Initializes the ScattData object from a given scatter and
     //!   multiplicity matrix.
@@ -112,10 +110,9 @@ class ScattData {
     //! @param gin Incoming energy group.
     //! @param gout Sampled outgoing energy group.
     //! @param i_gout Sampled outgoing energy group index.
-    //! @param prn_seeds Array of pseudorandom stream seeds
-    //! @param stream Pseudorandom stream index
+    //! @param prn_seed Pseudorandom number seed pointer
     void
-    sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seeds, int stream);
+    sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seed);
 
     //! \brief Provides a cross section value given certain parameters
     //!
@@ -166,7 +163,7 @@ class ScattDataLegendre: public ScattData {
     calc_f(int gin, int gout, double mu);
 
     void
-    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seeds, int stream);
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
 
     size_t
     get_order() {return dist[0][0].size() - 1;};
@@ -202,7 +199,7 @@ class ScattDataHistogram: public ScattData {
     calc_f(int gin, int gout, double mu);
 
     void
-    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seeds, int stream);
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
 
     size_t
     get_order() {return dist[0][0].size();};
@@ -243,7 +240,7 @@ class ScattDataTabular: public ScattData {
     calc_f(int gin, int gout, double mu);
 
     void
-    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seeds, int stream);
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
 
     size_t
     get_order() {return dist[0][0].size();};
diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
index 972df2f67..6b7c36796 100644
--- a/include/openmc/secondary_correlated.h
+++ b/include/openmc/secondary_correlated.h
@@ -38,10 +38,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-		  int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) 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 42c9130c1..b9f190635 100644
--- a/include/openmc/secondary_kalbach.h
+++ b/include/openmc/secondary_kalbach.h
@@ -29,10 +29,9 @@ 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
-  //! \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;
+  //! \param[inout] prn_seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) 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 0f8695c2e..b810d17d0 100644
--- a/include/openmc/secondary_nbody.h
+++ b/include/openmc/secondary_nbody.h
@@ -24,10 +24,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) 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 56cb16097..74bb3265f 100644
--- a/include/openmc/secondary_thermal.h
+++ b/include/openmc/secondary_thermal.h
@@ -31,10 +31,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) const override;
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
 };
@@ -54,10 +53,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) const override;
 private:
   double debye_waller_;
 };
@@ -78,10 +76,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) 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
@@ -103,10 +100,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) const override;
 private:
   const std::vector<double>& energy_; //!< Incident energies
   xt::xtensor<double, 2> energy_out_; //!< Outgoing energies for each incident energy
@@ -129,10 +125,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) const override;
 private:
   //! Secondary energy/angle distribution
   struct DistEnergySab {
diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h
index 0ac79527b..66ef834f1 100644
--- a/include/openmc/secondary_uncorrelated.h
+++ b/include/openmc/secondary_uncorrelated.h
@@ -29,10 +29,9 @@ 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
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
-  void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seeds,
-    int stream) const override;
+  //! \param[inout] prn_seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* prn_seed) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }
diff --git a/include/openmc/source.h b/include/openmc/source.h
index 6ea2b7ede..8a7c30a09 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -38,10 +38,9 @@ public:
   explicit SourceDistribution(pugi::xml_node node);
 
   //! Sample from the external source distribution
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
+  //! \param[inout] prn_seed Pseudorandom seed pointer
   //! \return Sampled site
-  Particle::Bank sample(uint64_t* prn_seed, int stream) const;
+  Particle::Bank sample(uint64_t* prn_seed) const;
 
   // Properties
   double strength() const { return strength_; }
@@ -62,9 +61,9 @@ extern "C" void initialize_source();
 
 //! Sample a site from all external source distributions in proportion to their
 //! source strength
-//! \param[inout] prn_seeds Array of pseudorandom seeds
+//! \param[inout] prn_seed Pseudorandom seed pointer
 //! \return Sampled source site
-Particle::Bank sample_external_source(uint64_t* prn_seeds);
+Particle::Bank sample_external_source(uint64_t* prn_seed);
 
 //! 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 3f17f4ec3..35fa84bd1 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -116,8 +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, uint64_t* prn_seeds,
-    int stream) const;
+  virtual Direction diffuse_reflect(Position r, Direction u,
+    uint64_t* prn_seed) const;
 
   //! Evaluate the equation describing the surface.
   //!
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index ed57e4a09..37ad8f4cb 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -64,10 +64,9 @@ 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
+  //! \param[inout] prn_seed Pseudorandom seed pointer
   void sample(const NuclideMicroXS& micro_xs, double E_in,
-              double* E_out, double* mu, uint64_t* prn_seeds, int stream);
+              double* E_out, double* mu, uint64_t* prn_seed);
 private:
   struct Reaction {
     // Default constructor
@@ -102,10 +101,9 @@ public:
   //! \param[out] i_temp corresponding temperature index
   //! \param[out] elastic Thermal elastic scattering cross section
   //! \param[out] inelastic Thermal inelastic scattering cross section
-  //! \param[inout] prn_seeds Array of pseudorandom seeds
-  //! \param[in] stream Pseudorandom stream index
+  //! \param[inout] prn_seed Pseudorandom seed pointer
   void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic,
-                    double* inelastic, uint64_t* prn_seeds, int stream) const;
+                    double* inelastic, uint64_t* prn_seed) const;
 
   //! Determine whether table applies to a particular nuclide
   //!
diff --git a/openmc/lib/math.py b/openmc/lib/math.py
index f1997b351..991df8f9c 100644
--- a/openmc/lib/math.py
+++ b/openmc/lib/math.py
@@ -5,6 +5,8 @@ from numpy.ctypeslib import ndpointer
 
 from . import _dll
 
+from random import randint
+
 
 _dll.t_percentile.restype = c_double
 _dll.t_percentile.argtypes = [c_double, c_int]
@@ -26,19 +28,19 @@ _dll.calc_zn_rad.argtypes = [c_int, c_double, ndpointer(c_double)]
 
 _dll.rotate_angle_c.restype = None
 _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
-                                POINTER(c_double), ndpointer(c_uint64), c_int]
+                                POINTER(c_double), POINTER(c_uint64)]
 _dll.maxwell_spectrum.restype = c_double
-_dll.maxwell_spectrum.argtypes = [c_double, ndpointer(c_uint64), c_int]
+_dll.maxwell_spectrum.argtypes = [c_double, POINTER(c_uint64)]
 
 _dll.watt_spectrum.restype = c_double
-_dll.watt_spectrum.argtypes = [c_double, c_double, ndpointer(c_uint64), c_int]
+_dll.watt_spectrum.argtypes = [c_double, c_double, POINTER(c_uint64)]
 
 _dll.broaden_wmp_polynomials.restype = None
 _dll.broaden_wmp_polynomials.argtypes = [c_double, c_double, c_int,
                                          ndpointer(c_double)]
 
 _dll.normal_variate.restype = c_double
-_dll.normal_variate.argtypes = [c_double, c_double, ndpointer(c_uint64), c_int]
+_dll.normal_variate.argtypes = [c_double, c_double, POINTER(c_uint64)]
 
 def t_percentile(p, df):
     """ Calculate the percentile of the Student's t distribution with a
@@ -185,7 +187,7 @@ def calc_zn_rad(n, rho):
     return zn_rad
 
 
-def rotate_angle(uvw0, mu, phi, prn_seeds, stream ):
+def rotate_angle(uvw0, mu, phi, prn_seed=None):
     """ Rotates direction cosines through a polar angle whose cosine is
     mu and through an azimuthal angle sampled uniformly.
 
@@ -197,10 +199,8 @@ def rotate_angle(uvw0, mu, phi, prn_seeds, stream ):
         Polar angle cosine to rotate
     phi : float
         Azimuthal angle; if None, one will be sampled uniformly
-    prn_seeds : iterable of int
-        PRNG seed array
-    stream : int
-        PRNG stream index
+    prn_seed : int
+        PRNG seed; if None, one will be generated randomly
 
     Returns
     -------
@@ -209,19 +209,21 @@ def rotate_angle(uvw0, mu, phi, prn_seeds, stream ):
 
     """
 
+    if prn_seed is None:
+        prn_seed = randint(0,10000000)
+
     uvw0_arr = np.array(uvw0, dtype=np.float64)
-    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
     if phi is None:
-        _dll.rotate_angle_c(uvw0_arr, mu, None, prn_seeds_arr, stream)
+        _dll.rotate_angle_c(uvw0_arr, mu, None, c_uint64(prn_seed))
     else:
-        _dll.rotate_angle_c(uvw0_arr, mu, c_double(phi), prn_seeds_arr, stream)
+        _dll.rotate_angle_c(uvw0_arr, mu, c_double(phi), c_uint64(prn_seed))
 
     uvw = uvw0_arr
 
     return uvw
 
 
-def maxwell_spectrum(T, prn_seeds, stream):
+def maxwell_spectrum(T, prn_seed=None):
     """ Samples an energy from the Maxwell fission distribution based
     on a direct sampling scheme.
 
@@ -229,10 +231,8 @@ def maxwell_spectrum(T, prn_seeds, stream):
     ----------
     T : float
         Spectrum parameter
-    prn_seeds : iterable of int
-        PRNG seed array
-    stream : int
-        PRNG stream index
+    prn_seed : int
+        PRNG seed; if None, one will be generated randomly
 
     Returns
     -------
@@ -240,13 +240,14 @@ def maxwell_spectrum(T, prn_seeds, stream):
         Sampled outgoing energy
 
     """
+    
+    if prn_seed is None:
+        prn_seed = randint(0,10000000)
 	
-    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
-
-    return _dll.maxwell_spectrum(T, prn_seeds_arr, stream)
+    return _dll.maxwell_spectrum(T, c_uint64(prn_seed))
 
 
-def watt_spectrum(a, b, prn_seeds, stream):
+def watt_spectrum(a, b, prn_seed=None):
     """ Samples an energy from the Watt energy-dependent fission spectrum.
 
     Parameters
@@ -255,10 +256,8 @@ def watt_spectrum(a, b, prn_seeds, stream):
         Spectrum parameter a
     b : float
         Spectrum parameter b
-    prn_seeds : iterable of int
-        PRNG seed array
-    stream : int
-        PRNG stream index
+    prn_seed : int
+        PRNG seed; if None, one will be generated randomly
 
     Returns
     -------
@@ -266,13 +265,14 @@ def watt_spectrum(a, b, prn_seeds, stream):
         Sampled outgoing energy
 
     """
-	
-    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
+    
+    if prn_seed is None:
+        prn_seed = randint(0,10000000)
 
-    return _dll.watt_spectrum(a, b, prn_seeds_arr, stream)
+    return _dll.watt_spectrum(a, b, c_uint64(prn_seed))
 
 
-def normal_variate(mean_value, std_dev, prn_seeds, stream):
+def normal_variate(mean_value, std_dev, prn_seed=None):
     """ Samples an energy from the Normal distribution.
 
     Parameters
@@ -281,10 +281,8 @@ def normal_variate(mean_value, std_dev, prn_seeds, stream):
         Mean of the Normal distribution
     std_dev : float
         Standard deviation of the normal distribution
-    prn_seeds : iterable of int
-        PRNG seed array
-    stream : int
-        PRNG stream index
+    prn_seed : int
+        PRNG seed; if None, one will be generated randomly
 
     Returns
     -------
@@ -292,10 +290,11 @@ def normal_variate(mean_value, std_dev, prn_seeds, stream):
         Sampled outgoing normally distributed value
 
     """
-	
-    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
+    
+    if prn_seed is None:
+        prn_seed = randint(0,10000000)
 
-    return _dll.normal_variate(mean_value, std_dev, prn_seeds_arr, stream)
+    return _dll.normal_variate(mean_value, std_dev, c_uint64(prn_seed))
 
 
 def broaden_wmp_polynomials(E, dopp, n):
diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp
index a99d118f0..85733eccc 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.cpp b/src/distribution.cpp
index 8083017b1..e48ce8baa 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -35,11 +35,11 @@ Discrete::Discrete(const double* x, const double* p, int n)
   normalize();
 }
 
-double Discrete::sample(uint64_t* prn_seeds, int stream) const
+double Discrete::sample(uint64_t* prn_seed) const
 {
   int n = x_.size();
   if (n > 1) {
-    double xi = prn(prn_seeds, stream);
+    double xi = prn(prn_seed);
     double c = 0.0;
     for (int i = 0; i < n; ++i) {
       c += p_[i];
@@ -74,9 +74,9 @@ Uniform::Uniform(pugi::xml_node node)
   b_ = params.at(1);
 }
 
-double Uniform::sample(uint64_t* prn_seeds, int stream) const
+double Uniform::sample(uint64_t* prn_seed) const
 {
-  return a_ + prn(prn_seeds, stream)*(b_ - a_);
+  return a_ + prn(prn_seed)*(b_ - a_);
 }
 
 //==============================================================================
@@ -88,9 +88,9 @@ Maxwell::Maxwell(pugi::xml_node node)
   theta_ = std::stod(get_node_value(node, "parameters"));
 }
 
-double Maxwell::sample(uint64_t* prn_seeds, int stream) const
+double Maxwell::sample(uint64_t* prn_seed) const
 {
-  return maxwell_spectrum(theta_, prn_seeds, stream);
+  return maxwell_spectrum(theta_, prn_seed);
 }
 
 //==============================================================================
@@ -108,9 +108,9 @@ Watt::Watt(pugi::xml_node node)
   b_ = params.at(1);
 }
 
-double Watt::sample(uint64_t* prn_seeds, int stream) const
+double Watt::sample(uint64_t* prn_seed) const
 {
-  return watt_spectrum(a_, b_, prn_seeds, stream);
+  return watt_spectrum(a_, b_, prn_seed);
 }
 
 //==============================================================================
@@ -127,9 +127,9 @@ Normal::Normal(pugi::xml_node node)
   std_dev_ = params.at(1);
 }
 
-double Normal::sample(uint64_t* prn_seeds, int stream) const
+double Normal::sample(uint64_t* prn_seed) const
 {
-  return normal_variate(mean_value_, std_dev_, prn_seeds, stream);
+  return normal_variate(mean_value_, std_dev_, prn_seed);
 }
 
 //==============================================================================
@@ -147,9 +147,9 @@ Muir::Muir(pugi::xml_node node)
   kt_ = params.at(2);
 }
 
-double Muir::sample(uint64_t* prn_seeds, int stream) const
+double Muir::sample(uint64_t* prn_seed) const
 {
-  return muir_spectrum(e0_, m_rat_, kt_, prn_seeds, stream);
+  return muir_spectrum(e0_, m_rat_, kt_, prn_seed);
 }
 
 //==============================================================================
@@ -220,10 +220,10 @@ void Tabular::init(const double* x, const double* p, std::size_t n, const double
   }
 }
 
-double Tabular::sample(uint64_t* prn_seeds, int stream) const
+double Tabular::sample(uint64_t* prn_seed) const
 {
   // Sample value of CDF
-  double c = prn(prn_seeds, stream);
+  double c = prn(prn_seed);
 
   // Find first CDF bin which is above the sampled value
   double c_i = c_[0];
@@ -263,11 +263,11 @@ double Tabular::sample(uint64_t* prn_seeds, int stream) const
 // Equiprobable implementation
 //==============================================================================
 
-double Equiprobable::sample(uint64_t* prn_seeds, int stream) const
+double Equiprobable::sample(uint64_t* prn_seed) const
 {
   std::size_t n = x_.size();
 
-  double r = prn(prn_seeds, stream);
+  double r = prn(prn_seed);
   int i = std::floor((n - 1)*r);
 
   double xl = x_[i];
diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index ec22bdd8c..076f82470 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -62,7 +62,7 @@ AngleDistribution::AngleDistribution(hid_t group)
   }
 }
 
-double AngleDistribution::sample(double E, uint64_t* prn_seeds, int stream) const
+double AngleDistribution::sample(double E, uint64_t* prn_seed) const
 {
   // Determine number of incoming energies
   auto n = energy_.size();
@@ -83,10 +83,10 @@ double AngleDistribution::sample(double E, uint64_t* prn_seeds, int stream) cons
   }
 
   // Sample between the ith and (i+1)th bin
-  if (r > prn(prn_seeds, stream)) ++i;
+  if (r > prn(prn_seed)) ++i;
 
   // Sample i-th distribution
-  double mu = distribution_[i]->sample(prn_seeds, stream);
+  double mu = distribution_[i]->sample(prn_seed);
 
   // Make sure mu is in range [-1,1] and return
   if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index cfb37abf3..3e1995145 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, uint64_t* prn_seeds, int stream) const
+double DiscretePhoton::sample(double E, uint64_t* prn_seed) 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, uint64_t* prn_seeds, int stream ) const
+double LevelInelastic::sample(double E, uint64_t* prn_seed ) const
 {
   return mass_ratio_*(E - threshold_);
 }
@@ -146,7 +146,7 @@ ContinuousTabular::ContinuousTabular(hid_t group)
   } // incoming energies
 }
 
-double ContinuousTabular::sample(double E, uint64_t* prn_seeds, int stream) const
+double ContinuousTabular::sample(double E, uint64_t* prn_seed) const
 {
   // Read number of interpolation regions and incoming energies
   bool histogram_interp;
@@ -177,7 +177,7 @@ double ContinuousTabular::sample(double E, uint64_t* prn_seeds, int stream) cons
   if (histogram_interp) {
     l = i;
   } else {
-    l = r > prn(prn_seeds, stream) ? i + 1 : i;
+    l = r > prn(prn_seed) ? i + 1 : i;
   }
 
   // Interpolation for energy E1 and EK
@@ -197,7 +197,7 @@ double ContinuousTabular::sample(double E, uint64_t* prn_seeds, int stream) cons
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
   n_discrete = distribution_[l].n_discrete;
-  double r1 = prn(prn_seeds, stream);
+  double r1 = prn(prn_seed);
   double c_k = distribution_[l].c[0];
   int k = 0;
   int end = n_energy_out - 2;
@@ -275,14 +275,14 @@ MaxwellEnergy::MaxwellEnergy(hid_t group)
   close_dataset(dset);
 }
 
-double MaxwellEnergy::sample(double E, uint64_t* prn_seeds, int stream) const
+double MaxwellEnergy::sample(double E, uint64_t* prn_seed) const
 {
   // Get temperature corresponding to incoming energy
   double theta = theta_(E);
 
   while (true) {
     // Sample maxwell fission spectrum
-    double E_out = maxwell_spectrum(theta, prn_seeds, stream);
+    double E_out = maxwell_spectrum(theta, prn_seed);
 
     // Accept energy based on restriction energy
     if (E_out <= E - u_) return E_out;
@@ -301,7 +301,7 @@ Evaporation::Evaporation(hid_t group)
   close_dataset(dset);
 }
 
-double Evaporation::sample(double E, uint64_t* prn_seeds, int stream) const
+double Evaporation::sample(double E, uint64_t* prn_seed) const
 {
   // Get temperature corresponding to incoming energy
   double theta = theta_(E);
@@ -313,7 +313,7 @@ double Evaporation::sample(double E, uint64_t* prn_seeds, int stream) const
   // density function
   double x;
   while (true) {
-    x = -std::log((1.0 - v*prn(prn_seeds, stream))*(1.0 - v*prn(prn_seeds, stream)));
+    x = -std::log((1.0 - v*prn(prn_seed))*(1.0 - v*prn(prn_seed)));
     if (x <= y) break;
   }
 
@@ -338,7 +338,7 @@ WattEnergy::WattEnergy(hid_t group)
   close_dataset(dset);
 }
 
-double WattEnergy::sample(double E, uint64_t* prn_seeds, int stream) const
+double WattEnergy::sample(double E, uint64_t* prn_seed) const
 {
   // Determine Watt parameters at incident energy
   double a = a_(E);
@@ -346,7 +346,7 @@ double WattEnergy::sample(double E, uint64_t* prn_seeds, int stream) const
 
   while (true) {
     // Sample energy-dependent Watt fission spectrum
-    double E_out = watt_spectrum(a, b, prn_seeds, stream);
+    double E_out = watt_spectrum(a, b, prn_seed);
 
     // Accept energy based on restriction energy
     if (E_out <= E - u_) return E_out;
diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp
index bbaa26de1..711234201 100644
--- a/src/distribution_multi.cpp
+++ b/src/distribution_multi.cpp
@@ -53,31 +53,31 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
   }
 }
 
-Direction PolarAzimuthal::sample(uint64_t* prn_seeds, int stream) const
+Direction PolarAzimuthal::sample(uint64_t* prn_seed) const
 {
   // Sample cosine of polar angle
-  double mu = mu_->sample(prn_seeds, stream);
+  double mu = mu_->sample(prn_seed);
   if (mu == 1.0) return u_ref_;
 
   // Sample azimuthal angle
-  double phi = phi_->sample(prn_seeds, stream);
+  double phi = phi_->sample(prn_seed);
 
   // If the reference direction is along the z-axis, rotate the aziumthal angle
   // to match spherical coordinate conventions.
   // 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, prn_seeds, stream);
+  return rotate_angle(u_ref_, mu, &phi, prn_seed);
 }
 
 //==============================================================================
 // Isotropic implementation
 //==============================================================================
 
-Direction Isotropic::sample(uint64_t* prn_seeds, int stream) const
+Direction Isotropic::sample(uint64_t* prn_seed) const
 {
-  double phi = 2.0*PI*prn(prn_seeds, stream);
-  double mu = 2.0*prn(prn_seeds, stream) - 1.0;
+  double phi = 2.0*PI*prn(prn_seed);
+  double mu = 2.0*prn(prn_seed) - 1.0;
   return {mu, std::sqrt(1.0 - mu*mu) * std::cos(phi),
       std::sqrt(1.0 - mu*mu) * std::sin(phi)};
 }
@@ -86,7 +86,7 @@ Direction Isotropic::sample(uint64_t* prn_seeds, int stream) const
 // Monodirectional implementation
 //==============================================================================
 
-Direction Monodirectional::sample(uint64_t* prn_seeds, int stream) const
+Direction Monodirectional::sample(uint64_t* prn_seed) const
 {
   return u_ref_;
 }
diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp
index e695322cd..ce2b7f34a 100644
--- a/src/distribution_spatial.cpp
+++ b/src/distribution_spatial.cpp
@@ -46,9 +46,9 @@ CartesianIndependent::CartesianIndependent(pugi::xml_node node)
   }
 }
 
-Position CartesianIndependent::sample(uint64_t* prn_seeds, int stream) const
+Position CartesianIndependent::sample(uint64_t* prn_seed) const
 {
-  return {x_->sample(prn_seeds, stream), y_->sample(prn_seeds, stream), z_->sample(prn_seeds, stream)};
+  return {x_->sample(prn_seed), y_->sample(prn_seed), z_->sample(prn_seed)};
 }
 
 //==============================================================================
@@ -107,11 +107,11 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
 
 }
 
-Position SphericalIndependent::sample(uint64_t* prn_seeds, int stream) const
+Position SphericalIndependent::sample(uint64_t* prn_seed) const
 {
-  double r = r_->sample(prn_seeds, stream);
-  double theta = theta_->sample(prn_seeds, stream);
-  double phi = phi_->sample(prn_seeds, stream);
+  double r = r_->sample(prn_seed);
+  double theta = theta_->sample(prn_seed);
+  double phi = phi_->sample(prn_seed);
   double x = r*sin(theta)*cos(phi) + origin_.x;
   double y =  r*sin(theta)*sin(phi) + origin_.y;
   double z = r*cos(theta) + origin_.z;
@@ -135,9 +135,9 @@ SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
   upper_right_ = Position{params[3], params[4], params[5]};
 }
 
-Position SpatialBox::sample(uint64_t* prn_seeds, int stream) const
+Position SpatialBox::sample(uint64_t* prn_seed) const
 {
-  Position xi {prn(prn_seeds, stream), prn(prn_seeds, stream), prn(prn_seeds, stream)};
+  Position xi {prn(prn_seed), prn(prn_seed), prn(prn_seed)};
   return lower_left_ + xi*(upper_right_ - lower_left_);
 }
 
@@ -157,7 +157,7 @@ SpatialPoint::SpatialPoint(pugi::xml_node node)
   r_ = Position{params.data()};
 }
 
-Position SpatialPoint::sample(uint64_t* prn_seeds, int stream) const
+Position SpatialPoint::sample(uint64_t* prn_seed) const
 {
   return r_;
 }
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index 4d4ac420e..cd0699a32 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -114,16 +114,16 @@ void synchronize_bank()
     fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
   }
 
-  // Create pseudorandom number streams
-  uint64_t prn_seeds[N_STREAMS];     // seeds
-  int stream = STREAM_TRACKING; // current RNG stream
+  // 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.
 
-  set_particle_seed(simulation::total_gen + overall_generation(), prn_seeds);
-  advance_prn_seed(start, prn_seeds, stream);
+  int64_t id = simulation::total_gen + overall_generation();
+  init_seed(id, &prn_seed, STREAM_TRACKING);
+  advance_prn_seed(start, &prn_seed);
 
   // Determine how many fission sites we need to sample from the source bank
   // and the probability for selecting a site.
@@ -158,7 +158,7 @@ void synchronize_bank()
     }
 
     // Randomly sample sites needed
-    if (prn(prn_seeds, stream) < p_sample) {
+    if (prn(&prn_seed) < p_sample) {
       temp_sites[index_temp] = site;
       ++index_temp;
     }
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 8db83609c..21afd696c 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -630,22 +630,22 @@ void calc_zn_rad(int n, double rho, double zn_rad[]) {
 }
 
 
-void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t* prn_seeds, int stream) {
-  Direction u = rotate_angle({uvw}, mu, phi, prn_seeds, stream);
+void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t* prn_seed) {
+  Direction u = rotate_angle({uvw}, mu, phi, prn_seed);
   uvw[0] = u.x;
   uvw[1] = u.y;
   uvw[2] = u.z;
 }
 
 
-Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* prn_seeds, int stream)
+Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* prn_seed)
 {
   // Sample azimuthal angle in [0,2pi) if none provided
   double phi_;
   if (phi != nullptr) {
     phi_ = (*phi);
   } else {
-    phi_ = 2.0*PI*prn(prn_seeds, stream);
+    phi_ = 2.0*PI*prn(prn_seed);
   }
 
   // Precompute factors to save flops
@@ -675,11 +675,11 @@ Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* prn_
 }
 
 
-double maxwell_spectrum(double T, uint64_t* prn_seeds, int stream) {
+double maxwell_spectrum(double T, uint64_t* prn_seed) {
   // Set the random numbers
-  double r1 = prn(prn_seeds, stream);
-  double r2 = prn(prn_seeds, stream);
-  double r3 = prn(prn_seeds, stream);
+  double r1 = prn(prn_seed);
+  double r2 = prn(prn_seed);
+  double r3 = prn(prn_seed);
 
   // determine cosine of pi/2*r
   double c = std::cos(PI / 2. * r3);
@@ -691,33 +691,33 @@ double maxwell_spectrum(double T, uint64_t* prn_seeds, int stream) {
 }
 
 
-double normal_variate(double mean, double standard_deviation, uint64_t* prn_seeds, int stream) {
+double normal_variate(double mean, double standard_deviation, uint64_t* prn_seed) {
   // perhaps there should be a limit to the number of resamples
   while ( true ) {
-    double v1 = 2 * prn(prn_seeds, stream) - 1.;
-    double v2 = 2 * prn(prn_seeds, stream) - 1.;
+    double v1 = 2 * prn(prn_seed) - 1.;
+    double v2 = 2 * prn(prn_seed) - 1.;
 
     double r = std::pow(v1, 2) + std::pow(v2, 2);
     double r2 = std::pow(r, 2);
     if (r2 < 1) {
       double z = std::sqrt(-2.0 * std::log(r2)/r2);
-      z *= (prn(prn_seeds, stream) <= 0.5) ? v1 : v2;
+      z *= (prn(prn_seed) <= 0.5) ? v1 : v2;
       return mean + standard_deviation*z;
     }
   }
 }
 
-double muir_spectrum(double e0, double m_rat, double kt, uint64_t* prn_seeds, int stream) {
+double muir_spectrum(double e0, double m_rat, double kt, uint64_t* prn_seed) {
   // note sigma here is a factor of 2 shy of equation
   // 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
   double sigma = std::sqrt(2.*e0*kt/m_rat);
-  return normal_variate(e0, sigma, prn_seeds, stream);
+  return normal_variate(e0, sigma, prn_seed);
 }
 
 
-double watt_spectrum(double a, double b, uint64_t* prn_seeds, int stream) {
-  double w = maxwell_spectrum(a, prn_seeds, stream);
-  double E_out = w + 0.25 * a * a * b + (2. * prn(prn_seeds, stream) - 1.) * std::sqrt(a * a * b * w);
+double watt_spectrum(double a, double b, uint64_t* prn_seed) {
+  double w = maxwell_spectrum(a, prn_seed);
+  double E_out = w + 0.25 * a * a * b + (2. * prn(prn_seed) - 1.) * std::sqrt(a * a * b * w);
 
   return E_out;
 }
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index 9ead9b807..e2195d512 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -529,7 +529,7 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
 //==============================================================================
 
 void
-Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seeds, int stream)
+Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed)
 {
   // This method assumes that the temperature and angle indices are set
 #ifdef _OPENMP
@@ -544,8 +544,8 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seeds, in
   double prob_prompt = xs_t->prompt_nu_fission(cache[tid].a, gin);
 
   // sample random numbers
-  double xi_pd = prn(prn_seeds, stream) * nu_fission;
-  double xi_gout = prn(prn_seeds, stream);
+  double xi_pd = prn(prn_seed) * nu_fission;
+  double xi_gout = prn(prn_seed);
 
   // Select whether the neutron is prompt or delayed
   if (xi_pd <= prob_prompt) {
@@ -585,8 +585,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seeds, in
 //==============================================================================
 
 void
-Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seeds,
-  int stream)
+Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed)
 {
   // This method assumes that the temperature and angle indices are set
   // Sample the data
@@ -595,7 +594,7 @@ Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_
 #else
   int tid = 0;
 #endif
-  xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, prn_seeds, stream);
+  xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, prn_seed);
 }
 
 //==============================================================================
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index b5da3f60a..bc0fed2e0 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);
@@ -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.prn_seeds_[p.stream_]);
   p.stream_ = STREAM_TRACKING;
 
   int i_low = 0;
diff --git a/src/particle.cpp b/src/particle.cpp
index 98ee8920b..d3db2e84b 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(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 17eb5a3e6..b7d94dd80 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_);
+  init_particle_seeds(particle_seed, p.prn_seeds_);
 
   // Transport neutron
   p.transport();
diff --git a/src/photon.cpp b/src/photon.cpp
index 7db4c3f29..6cf43cbc5 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -290,12 +290,12 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
 }
 
 void PhotonInteraction::compton_scatter(double alpha, bool doppler,
-  double* alpha_out, double* mu, int* i_shell, uint64_t* prn_seeds, int stream) const
+  double* alpha_out, double* mu, int* i_shell, uint64_t* prn_seed) const
 {
   double form_factor_xmax = 0.0;
   while (true) {
     // Sample Klein-Nishina distribution for trial energy and angle
-    std::tie(*alpha_out, *mu) = klein_nishina(alpha, prn_seeds, stream);
+    std::tie(*alpha_out, *mu) = klein_nishina(alpha, prn_seed);
 
     // Note that the parameter used here does not correspond exactly to the
     // momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the
@@ -309,10 +309,10 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler,
     }
 
     // Perform rejection on form factor
-    if (prn(prn_seeds, stream) < form_factor_x / form_factor_xmax) {
+    if (prn(prn_seed) < form_factor_x / form_factor_xmax) {
       if (doppler) {
         double E_out;
-        this->compton_doppler(alpha, *mu, &E_out, i_shell, prn_seeds, stream);
+        this->compton_doppler(alpha, *mu, &E_out, i_shell, prn_seed);
         *alpha_out = E_out/MASS_ELECTRON_EV;
       } else {
         *i_shell = -1;
@@ -323,14 +323,14 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler,
 }
 
 void PhotonInteraction::compton_doppler(double alpha, double mu,
-  double* E_out, int* i_shell, uint64_t* prn_seeds, int stream) const
+  double* E_out, int* i_shell, uint64_t* prn_seed) const
 {
   auto n = data::compton_profile_pz.size();
 
   int shell; // index for shell
   while (true) {
     // Sample electron shell
-    double rn = prn(prn_seeds, stream);
+    double rn = prn(prn_seed);
     double c = 0.0;
     for (shell = 0; shell < electron_pdf_.size(); ++shell) {
       c += electron_pdf_(shell);
@@ -377,7 +377,7 @@ void PhotonInteraction::compton_doppler(double alpha, double mu,
     }
 
     // Sample value on bounded cdf
-    c = prn(prn_seeds, stream)*c_max;
+    c = prn(prn_seed)*c_max;
 
     // Determine pz corresponding to sampled cdf value
     auto cdf_shell = xt::view(profile_cdf_, shell, xt::all());
@@ -418,7 +418,7 @@ void PhotonInteraction::compton_doppler(double alpha, double mu,
     if (E_out1 > 0.0) {
       if (E_out2 > 0.0) {
         // If both are positive, pick one at random
-        *E_out = prn(prn_seeds, stream) < 0.5 ? E_out1 : E_out2;
+        *E_out = prn(prn_seed) < 0.5 ? E_out1 : E_out2;
       } else {
         *E_out = E_out1;
       }
@@ -496,7 +496,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const
   xs.last_E = p.E_;
 }
 
-double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seeds, int stream) const
+double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seed) const
 {
   double mu;
   while (true) {
@@ -507,7 +507,7 @@ double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seeds, in
     double F_max = coherent_int_form_factor_(x2_max);
 
     // Sample cumulative distribution
-    double F = prn(prn_seeds, stream)*F_max;
+    double F = prn(prn_seed)*F_max;
 
     // Determine x^2 corresponding to F
     const auto& x {coherent_int_form_factor_.x()};
@@ -519,14 +519,14 @@ double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seeds, in
     // Calculate mu
     mu = 1.0 - 2.0*x2/x2_max;
 
-    if (prn(prn_seeds, stream) < 0.5*(1.0 + mu*mu)) break;
+    if (prn(prn_seed) < 0.5*(1.0 + mu*mu)) break;
   }
   return mu;
 }
 
 void PhotonInteraction::pair_production(double alpha, double* E_electron,
-  double* E_positron, double* mu_electron, double* mu_positron, uint64_t* prn_seeds,
-  int stream) const
+  double* E_positron, double* mu_electron, double* mu_positron,
+  uint64_t* prn_seed) const
 {
   constexpr double r[] {
     122.81, 73.167, 69.228, 67.301, 64.696, 61.228,
@@ -593,12 +593,12 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron,
   double u2 = phi2_max;
   double e;
   while (true) {
-    double rn = prn(prn_seeds, stream);
+    double rn = prn(prn_seed);
 
     // Sample the index i in (1, 2) using the point probabilities
     // p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2)
     int i;
-    if (prn(prn_seeds, stream) < u1/(u1 + u2)) {
+    if (prn(prn_seed) < u1/(u1 + u2)) {
       i = 1;
 
       // Sample e from pi_1 using the inverse transform method
@@ -619,10 +619,10 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron,
     t3 = b*b*(4.0 - 4.0*t2 - 3.0*std::log(1.0 + 1.0/(b*b)));
     if (i == 1) {
       double phi1 = 7.0/3.0 - t1 - 6.0*t2 - t3 + t4;
-      if (prn(prn_seeds, stream) <= phi1/phi1_max) break;
+      if (prn(prn_seed) <= phi1/phi1_max) break;
     } else {
       double phi2 = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4;
-      if (prn(prn_seeds, stream) <= phi2/phi2_max) break;
+      if (prn(prn_seed) <= phi2/phi2_max) break;
     }
   }
 
@@ -635,13 +635,13 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron,
   // p(mu) = C/(1 - beta*mu)^2 using the inverse transform method.
   double beta = std::sqrt(*E_electron*(*E_electron + 2.0*MASS_ELECTRON_EV))
     / (*E_electron + MASS_ELECTRON_EV)  ;
-  double rn = 2.0*prn(prn_seeds, stream) - 1.0;
+  double rn = 2.0*prn(prn_seed) - 1.0;
   *mu_electron = (rn + beta)/(rn*beta + 1.0);
 
   // Sample the scattering angle of the positron
   beta = std::sqrt(*E_positron*(*E_positron + 2.0*MASS_ELECTRON_EV))
     / (*E_positron + MASS_ELECTRON_EV);
-  rn = 2.0*prn(prn_seeds, stream) - 1.0;
+  rn = 2.0*prn(prn_seed) - 1.0;
   *mu_positron = (rn + beta)/(rn*beta + 1.0);
 }
 
@@ -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);
@@ -711,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
 // Non-member functions
 //==============================================================================
 
-std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seeds, int stream)
+std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seed)
 {
   double alpha_out, mu;
   double beta = 1.0 + 2.0*alpha;
@@ -720,19 +720,19 @@ std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seeds, int s
     double t = beta/(beta + 8.0);
     double x;
     while (true) {
-      if (prn(prn_seeds, stream) < t) {
+      if (prn(prn_seed) < t) {
         // Left branch of flow chart
-        double r = 2.0*prn(prn_seeds, stream);
+        double r = 2.0*prn(prn_seed);
         x = 1.0 + alpha*r;
-        if (prn(prn_seeds, stream) < 4.0/x*(1.0 - 1.0/x)) {
+        if (prn(prn_seed) < 4.0/x*(1.0 - 1.0/x)) {
           mu = 1 - r;
           break;
         }
       } else {
         // Right branch of flow chart
-        x = beta/(1.0 + 2.0*alpha*prn(prn_seeds, stream));
+        x = beta/(1.0 + 2.0*alpha*prn(prn_seed));
         mu = 1.0 + (1.0 - x)/alpha;
-        if (prn(prn_seeds, stream) < 0.5*(mu*mu + 1.0/x)) break;
+        if (prn(prn_seed) < 0.5*(mu*mu + 1.0/x)) break;
       }
     }
     alpha_out = alpha/x;
@@ -740,24 +740,24 @@ std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seeds, int s
   } else {
     // Koblinger's direct method
     double gamma = 1.0 - std::pow(beta, -2);
-    double s = prn(prn_seeds, stream)*(4.0/alpha + 0.5*gamma +
+    double s = prn(prn_seed)*(4.0/alpha + 0.5*gamma +
       (1.0 - (1.0 + beta)/(alpha*alpha))*std::log(beta));
     if (s <= 2.0/alpha) {
       // For first term, x = 1 + 2ar
       // Therefore, a' = a/(1 + 2ar)
-      alpha_out = alpha/(1.0 + 2.0*alpha*prn(prn_seeds, stream));
+      alpha_out = alpha/(1.0 + 2.0*alpha*prn(prn_seed));
     } else if (s <= 4.0/alpha) {
       // For third term, x = beta/(1 + 2ar)
       // Therefore, a' = a(1 + 2ar)/beta
-      alpha_out = alpha*(1.0 + 2.0*alpha*prn(prn_seeds, stream))/beta;
+      alpha_out = alpha*(1.0 + 2.0*alpha*prn(prn_seed))/beta;
     } else if (s <= 4.0/alpha + 0.5*gamma) {
       // For fourth term, x = 1/sqrt(1 - gamma*r)
       // Therefore, a' = a*sqrt(1 - gamma*r)
-      alpha_out = alpha*std::sqrt(1.0 - gamma*prn(prn_seeds, stream));
+      alpha_out = alpha*std::sqrt(1.0 - gamma*prn(prn_seed));
     } else {
       // For third term, x = beta^r
       // Therefore, a' = a/beta^r
-      alpha_out = alpha/std::pow(beta, prn(prn_seeds, stream));
+      alpha_out = alpha/std::pow(beta, prn(prn_seed));
     }
 
     // Calculate cosine of scattering angle based on basic relation
diff --git a/src/physics.cpp b/src/physics.cpp
index 84b318f84..53e4a5ac0 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->prn_seeds_ + p->stream_);
     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->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,15 +742,15 @@ 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,
-  Direction v_neut, double xs_eff, double kT, uint64_t* prn_seeds, int stream)
+  Direction v_neut, double xs_eff, double kT, uint64_t* prn_seed)
 {
   // check if nuclide is a resonant scatterer
   ResScatMethod sampling_method;
@@ -782,7 +782,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
   case ResScatMethod::cxs:
 
     // sample target velocity with the constant cross section (cxs) approx.
-    return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream);
+    return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed);
 
   case ResScatMethod::dbrc:
   case ResScatMethod::rvs: {
@@ -816,7 +816,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
     if (i_E_up == i_E_low) {
       // Handle degenerate case -- if the upper/lower bounds occur for the same
       // index, then using cxs is probably a good approximation
-      return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream);
+      return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed);
     }
 
     if (sampling_method == ResScatMethod::dbrc) {
@@ -840,7 +840,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
         Direction v_target;
         while (true) {
           // sample target velocity with the constant cross section (cxs) approx.
-          v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream);
+          v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed);
           Direction v_rel = v_neut - v_target;
           E_rel = v_rel.dot(v_rel);
           if (E_rel < E_up) break;
@@ -849,7 +849,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
         // perform Doppler broadening rejection correction (dbrc)
         double xs_0K = nuc->elastic_xs_0K(E_rel);
         double R = xs_0K / xs_max;
-        if (prn(prn_seeds, stream) < R) return v_target;
+        if (prn(prn_seed) < R) return v_target;
       }
 
     } else if (sampling_method == ResScatMethod::rvs) {
@@ -869,10 +869,10 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
 
       while (true) {
         // directly sample Maxwellian
-        double E_t = -kT * std::log(prn(prn_seeds, stream));
+        double E_t = -kT * std::log(prn(prn_seed));
 
         // sample a relative energy using the xs cdf
-        double cdf_rel = cdf_low + prn(prn_seeds, stream)*(cdf_up - cdf_low);
+        double cdf_rel = cdf_low + prn(prn_seed)*(cdf_up - cdf_low);
         int i_E_rel = lower_bound_index(&nuc->xs_cdf_[i_E_low-1],
           &nuc->xs_cdf_[i_E_up+1], cdf_rel);
         double E_rel = nuc->energy_0K_[i_E_low + i_E_rel];
@@ -890,7 +890,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
         if (std::abs(mu) < 1.0) {
           // set and accept target velocity
           E_t /= nuc->awr_;
-          return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, prn_seeds, stream);
+          return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, prn_seed);
         }
       }
     }
@@ -901,7 +901,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
 }
 
 Direction
-sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_t* prn_seeds, int stream)
+sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_t* prn_seed)
 {
   double beta_vn = std::sqrt(awr * E / kT);
   double alpha = 1.0/(1.0 + std::sqrt(PI)*beta_vn/2.0);
@@ -910,10 +910,10 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
   double mu;
   while (true) {
     // Sample two random numbers
-    double r1 = prn(prn_seeds, stream);
-    double r2 = prn(prn_seeds, stream);
+    double r1 = prn(prn_seed);
+    double r2 = prn(prn_seed);
 
-    if (prn(prn_seeds, stream) < alpha) {
+    if (prn(prn_seed) < alpha) {
       // With probability alpha, we sample the distribution p(y) =
       // y*e^(-y). This can be done with sampling scheme C45 frmo the Monte
       // Carlo sampler
@@ -925,7 +925,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
       // e^(-y^2). This can be done with sampling scheme C61 from the Monte
       // Carlo sampler
 
-      double c = std::cos(PI/2.0 * prn(prn_seeds, stream));
+      double c = std::cos(PI/2.0 * prn(prn_seed));
       beta_vt_sq = -std::log(r1) - std::log(r2)*c*c;
     }
 
@@ -933,14 +933,14 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
     double beta_vt = std::sqrt(beta_vt_sq);
 
     // Sample cosine of angle between neutron and target velocity
-    mu = 2.0*prn(prn_seeds, stream) - 1.0;
+    mu = 2.0*prn(prn_seed) - 1.0;
 
     // Determine rejection probability
     double accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq -
       2*beta_vn*beta_vt*mu) / (beta_vn + beta_vt);
 
     // Perform rejection sampling on vt and mu
-    if (prn(prn_seeds, stream) < accept_prob) break;
+    if (prn(prn_seed) < accept_prob) break;
   }
 
   // Determine speed of target nucleus
@@ -948,19 +948,19 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
 
   // Determine velocity vector of target nucleus based on neutron's velocity
   // and the sampled angle between them
-  return vt * rotate_angle(u, mu, nullptr, prn_seeds, stream);
+  return vt * rotate_angle(u, mu, nullptr, prn_seed);
 }
 
-void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, uint64_t* prn_seeds, int stream)
+void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, uint64_t* prn_seed)
 {
   // Sample cosine of angle -- fission neutrons are always emitted
   // isotropically. Sometimes in ACE data, fission reactions actually have
   // an angular distribution listed, but for those that do, it's simply just
   // a uniform distribution in mu
-  double mu = 2.0 * prn(prn_seeds, stream) - 1.0;
+  double mu = 2.0 * prn(prn_seed) - 1.0;
 
   // Sample azimuthal angle uniformly in [0,2*pi)
-  double phi = 2.0*PI*prn(prn_seeds, stream);
+  double phi = 2.0*PI*prn(prn_seed);
   site->u.x = mu;
   site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi);
   site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi);
@@ -971,12 +971,12 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part
   double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed);
   double beta = nu_d / nu_t;
 
-  if (prn(prn_seeds, stream) < beta) {
+  if (prn(prn_seed) < beta) {
     // ====================================================================
     // DELAYED NEUTRON SAMPLED
 
     // sampled delayed precursor group
-    double xi = prn(prn_seeds, stream)*nu_d;
+    double xi = prn(prn_seed)*nu_d;
     double prob = 0.0;
     int group;
     for (group = 1; group < nuc->n_precursor_; ++group) {
@@ -1000,7 +1000,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part
     while (true) {
       // sample from energy/angle distribution -- note that mu has already been
       // sampled above and doesn't need to be resampled
-      rx->products_[group].sample(E_in, site->E, mu, prn_seeds, stream);
+      rx->products_[group].sample(E_in, site->E, mu, prn_seed);
 
       // resample if energy is greater than maximum neutron energy
       constexpr int neutron = static_cast<int>(Particle::Type::neutron);
@@ -1025,7 +1025,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part
     // sample from prompt neutron energy distribution
     int n_sample = 0;
     while (true) {
-      rx->products_[0].sample(E_in, site->E, mu, prn_seeds, stream);
+      rx->products_[0].sample(E_in, site->E, mu, prn_seed);
 
       // resample if energy is greater than maximum neutron energy
       constexpr int neutron = static_cast<int>(Particle::Type::neutron);
@@ -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 e5b43e223..c26ebdfd2 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 6898e1e0e..3313f1034 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/plot.cpp b/src/plot.cpp
index 07c612952..dbc1ccb9e 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -80,8 +80,7 @@ namespace model {
 
 std::vector<Plot> plots;
 std::unordered_map<int, int> plot_map;
-uint64_t plotter_prn_seeds[N_STREAMS] = {1, 2, 3, 4, 5, 6};
-int plotter_stream = STREAM_TRACKING;
+uint64_t plotter_prn_seed = 1;
 
 } // namespace model
 
@@ -961,9 +960,9 @@ voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
 }
 
 RGBColor random_color(void) {
-  return {int(prn(model::plotter_prn_seeds, model::plotter_stream)*255),
-          int(prn(model::plotter_prn_seeds, model::plotter_stream)*255),
-          int(prn(model::plotter_prn_seeds, model::plotter_stream)*255)};
+  return {int(prn(&model::plotter_prn_seed)*255),
+          int(prn(&model::plotter_prn_seed)*255),
+          int(prn(&model::plotter_prn_seed)*255)};
 }
 
 extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp
index 6483d86f5..2b5f66d59 100644
--- a/src/random_lcg.cpp
+++ b/src/random_lcg.cpp
@@ -32,15 +32,15 @@ constexpr double   prn_norm   {1.0 / prn_mod};           // 2^-63
 //==============================================================================
 
 extern "C" double
-prn(uint64_t* prn_seeds, int stream)
+prn(uint64_t* prn_seed)
 {
   // This algorithm uses bit-masking to find the next integer(8) value to be
   // used to calculate the random number.
-  prn_seeds[stream] = (prn_mult*prn_seeds[stream] + prn_add) & prn_mask;
+  *prn_seed = (prn_mult * (*prn_seed) + prn_add) & prn_mask;
 
   // Once the integer is calculated, we just need to divide by 2**m,
   // represented here as multiplying by a pre-calculated factor
-  return prn_seeds[stream] * prn_norm;
+  return (*prn_seed) * prn_norm;
 }
 
 //==============================================================================
@@ -48,17 +48,27 @@ prn(uint64_t* prn_seeds, int stream)
 //==============================================================================
 
 extern "C" double
-future_prn(int64_t n, uint64_t* prn_seeds, int stream)
+future_prn(int64_t n, uint64_t prn_seed)
 {
-  return future_seed(static_cast<uint64_t>(n), prn_seeds[stream]) * prn_norm;
+  return future_seed(static_cast<uint64_t>(n), prn_seed) * prn_norm;
 }
 
 //==============================================================================
-// SET_PARTICLE_SEED
+// SET_SEED
 //==============================================================================
 
 extern "C" void
-set_particle_seed(int64_t id, uint64_t* prn_seeds)
+init_seed(int64_t id, uint64_t* prn_seed, int offset)
+{
+  *prn_seed = future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + offset);
+}
+
+//==============================================================================
+// SET_PARTICLE_SEEDS
+//==============================================================================
+
+extern "C" void
+init_particle_seeds(int64_t id, uint64_t* prn_seeds)
 {
   for (int i = 0; i < N_STREAMS; i++) {
     prn_seeds[i] = future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + i);
@@ -70,9 +80,9 @@ set_particle_seed(int64_t id, uint64_t* prn_seeds)
 //==============================================================================
 
 extern "C" void
-advance_prn_seed(int64_t n, uint64_t* prn_seeds, int stream)
+advance_prn_seed(int64_t n, uint64_t* prn_seed)
 {
-  prn_seeds[stream] = future_seed(static_cast<uint64_t>(n), prn_seeds[stream]);
+  *prn_seed = future_seed(static_cast<uint64_t>(n), *prn_seed);
 }
 
 //==============================================================================
@@ -80,7 +90,7 @@ advance_prn_seed(int64_t n, uint64_t* prn_seeds, int stream)
 //==============================================================================
 
 uint64_t
-future_seed(uint64_t n, uint64_t seed)
+future_seed(uint64_t n, uint64_t prn_seed)
 {
   // Make sure nskip is less than 2^M.
   n &= prn_mask;
@@ -111,7 +121,7 @@ future_seed(uint64_t n, uint64_t seed)
   }
 
   // With G and C, we can now find the new seed.
-  return (g_new * seed + c_new) & prn_mask;
+  return (g_new * prn_seed + c_new) & prn_mask;
 }
 
 //==============================================================================
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
index fd2709b2d..d92d14517 100644
--- a/src/reaction_product.cpp
+++ b/src/reaction_product.cpp
@@ -77,25 +77,25 @@ ReactionProduct::ReactionProduct(hid_t group)
 }
 
 void ReactionProduct::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) const
 {
   auto n = applicability_.size();
   if (n > 1) {
     double prob = 0.0;
-    double c = prn(prn_seeds, stream);
+    double c = prn(prn_seed);
     for (int i = 0; i < n; ++i) {
       // Determine probability that i-th energy distribution is sampled
       prob += applicability_[i](E_in);
 
       // If i-th distribution is sampled, sample energy from the distribution
       if (c <= prob) {
-        distribution_[i]->sample(E_in, E_out, mu, prn_seeds, stream);
+        distribution_[i]->sample(E_in, E_out, mu, prn_seed);
         break;
       }
     }
   } else {
     // If only one distribution is present, go ahead and sample it
-    distribution_[0]->sample(E_in, E_out, mu, prn_seeds, stream);
+    distribution_[0]->sample(E_in, E_out, mu, prn_seed);
   }
 }
 
diff --git a/src/scattdata.cpp b/src/scattdata.cpp
index 89509fabe..5d06375e5 100644
--- a/src/scattdata.cpp
+++ b/src/scattdata.cpp
@@ -167,10 +167,10 @@ ScattData::base_combine(size_t max_order,
 //==============================================================================
 
 void
-ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seeds, int stream)
+ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seed)
 {
   // Sample the outgoing group
-  double xi = prn(prn_seeds, stream);
+  double xi = prn(prn_seed);
   double prob = 0.;
   i_gout = 0;
   for (gout = gmin[gin]; gout < gmax[gin]; ++gout) {
@@ -348,21 +348,21 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu)
 
 void
 ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt,
-                          uint64_t* prn_seeds, int stream)
+                          uint64_t* prn_seed)
 {
   // Sample the outgoing energy using the base-class method
   int i_gout;
-  sample_energy(gin, gout, i_gout, prn_seeds, stream);
+  sample_energy(gin, gout, i_gout, prn_seed);
 
   // Now we can sample mu using the scattering kernel using rejection
   // sampling from a rectangular bounding box
   double M = max_val[gin][i_gout];
   int samples;
   for (samples = 0; samples < MAX_SAMPLE; ++samples) {
-    mu = 2. * prn(prn_seeds, stream) - 1.;
+    mu = 2. * prn(prn_seed) - 1.;
     double f = calc_f(gin, gout, mu);
     if (f > 0.) {
-      double u = prn(prn_seeds, stream) * M;
+      double u = prn(prn_seed) * M;
       if (u <= f) break;
     }
   }
@@ -537,14 +537,14 @@ ScattDataHistogram::calc_f(int gin, int gout, double mu)
 
 void
 ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt,
-                           uint64_t* prn_seeds, int stream)
+                           uint64_t* prn_seed)
 {
   // Sample the outgoing energy using the base-class method
   int i_gout;
-  sample_energy(gin, gout, i_gout, prn_seeds, stream);
+  sample_energy(gin, gout, i_gout, prn_seed);
 
   // Determine the outgoing cosine bin
-  double xi = prn(prn_seeds, stream);
+  double xi = prn(prn_seed);
 
   int imu;
   if (xi < dist[gin][i_gout][0]) {
@@ -556,7 +556,7 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt,
   }
 
   // Randomly select mu within the imu bin
-  mu = prn(prn_seeds, stream) * dmu + this->mu[imu];
+  mu = prn(prn_seed) * dmu + this->mu[imu];
 
   if (mu < -1.) {
     mu = -1.;
@@ -741,15 +741,15 @@ ScattDataTabular::calc_f(int gin, int gout, double mu)
 
 void
 ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt,
-                         uint64_t* prn_seeds, int stream)
+                         uint64_t* prn_seed)
 {
   // Sample the outgoing energy using the base-class method
   int i_gout;
-  sample_energy(gin, gout, i_gout, prn_seeds, stream);
+  sample_energy(gin, gout, i_gout, prn_seed);
 
   // Determine the outgoing cosine bin
   int NP = this->mu.shape()[0];
-  double xi = prn(prn_seeds, stream);
+  double xi = prn(prn_seed);
 
   double c_k = dist[gin][i_gout][0];
   int k;
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
index 3973c7596..81e1ec1f9 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -153,14 +153,14 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
 }
 
 void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) 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(prn_seeds, stream) - 1.0;
+  mu = 2.0*prn(prn_seed) - 1.0;
   // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
   // Find energy bin and calculate interpolation factor -- if the energy is
@@ -180,7 +180,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
   }
 
   // Sample between the ith and [i+1]th bin
-  int l = r > prn(prn_seeds, stream) ? i + 1 : i;
+  int l = r > prn(prn_seed) ? i + 1 : i;
 
   // Interpolation for energy E1 and EK
   int n_energy_out = distribution_[i].e_out.size();
@@ -199,7 +199,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
   n_discrete = distribution_[l].n_discrete;
-  double r1 = prn(prn_seeds, stream);
+  double r1 = prn(prn_seed);
   double c_k = distribution_[l].c[0];
   int k = 0;
   int end = n_energy_out - 2;
@@ -260,9 +260,9 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
 
   // Find correlated angular distribution for closest outgoing energy bin
   if (r1 - c_k < c_k1 - r1) {
-    mu = distribution_[l].angle[k]->sample(prn_seeds, stream);
+    mu = distribution_[l].angle[k]->sample(prn_seed);
   } else {
-    mu = distribution_[l].angle[k + 1]->sample(prn_seeds, stream);
+    mu = distribution_[l].angle[k + 1]->sample(prn_seed);
   }
 }
 
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index 628b1392f..b0f9854ac 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, uint64_t* prn_seeds, int stream) const
+void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* prn_seed) 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(prn_seeds, stream) - 1.0;
+  mu = 2.0*prn(prn_seed) - 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, uint64_t* prn_s
   }
 
   // Sample between the ith and [i+1]th bin
-  int l = r > prn(prn_seeds, stream) ? i + 1 : i;
+  int l = r > prn(prn_seed) ? 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, uint64_t* prn_s
   // Determine outgoing energy bin
   n_energy_out = distribution_[l].e_out.size();
   n_discrete = distribution_[l].n_discrete;
-  double r1 = prn(prn_seeds, stream);
+  double r1 = prn(prn_seed);
   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, uint64_t* prn_s
   }
 
   // Sampled correlated angle from Kalbach-Mann parameters
-  if (prn(prn_seeds, stream) > km_r) {
-    double T = (2.0*prn(prn_seeds, stream) - 1.0) * std::sinh(km_a);
+  if (prn(prn_seed) > km_r) {
+    double T = (2.0*prn(prn_seed) - 1.0) * std::sinh(km_a);
     mu = std::log(T + std::sqrt(T*T + 1.0))/km_a;
   } else {
-    double r1 = prn(prn_seeds, stream);
+    double r1 = prn(prn_seed);
     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 a76230cc0..46e56145b 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -22,38 +22,38 @@ NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)
 }
 
 void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) const
 {
   // By definition, the distribution of the angle is isotropic for an N-body
   // phase space distribution
-  mu = 2.0*prn(prn_seeds, stream) - 1.0;
+  mu = 2.0*prn(prn_seed) - 1.0;
 
   // Determine E_max parameter
   double Ap = mass_ratio_;
   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, prn_seeds, stream);
+  double x = maxwell_spectrum(1.0, prn_seed);
 
   double y;
   double r1, r2, r3, r4, r5, r6;
   switch (n_bodies_) {
   case 3:
-    y = maxwell_spectrum(1.0, prn_seeds, stream);
+    y = maxwell_spectrum(1.0, prn_seed);
     break;
   case 4:
-    r1 = prn(prn_seeds, stream);
-    r2 = prn(prn_seeds, stream);
-    r3 = prn(prn_seeds, stream);
+    r1 = prn(prn_seed);
+    r2 = prn(prn_seed);
+    r3 = prn(prn_seed);
     y = -std::log(r1*r2*r3);
     break;
   case 5:
-    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);
+    r1 = prn(prn_seed);
+    r2 = prn(prn_seed);
+    r3 = prn(prn_seed);
+    r4 = prn(prn_seed);
+    r5 = prn(prn_seed);
+    r6 = prn(prn_seed);
     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 a5190d0b2..13b5746fc 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -32,8 +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_seed) const
 {
   // Get index and interpolation factor for elastic grid
   int i;
@@ -43,7 +43,7 @@ CoherentElasticAE::sample(double E_in, double& E_out, double& mu, uint64_t* prn_
 
   // Sample a Bragg edge between 1 and i
   const auto& factors = xs_.factors();
-  double prob = prn(prn_seeds, stream) * factors[i+1];
+  double prob = prn(prn_seed) * factors[i+1];
   int k = 0;
   if (prob >= factors.front()) {
     k = lower_bound_index(factors.begin(), factors.begin() + (i+1), prob);
@@ -67,11 +67,11 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group)
 
 void
 IncoherentElasticAE::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) 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(prn_seeds, stream)*(std::exp(2.0*c) - 1))/c - 1.0;
+  mu = std::log(1.0 + prn(prn_seed)*(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;
@@ -90,7 +90,7 @@ IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group,
 
 void
 IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) const
 {
   // Get index and interpolation factor for elastic grid
   int i;
@@ -101,7 +101,7 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
   // incoming energies.
 
   // Sample outgoing cosine bin
-  int k = prn(prn_seeds, stream) * mu_out_.shape()[1];
+  int k = prn(prn_seed) * mu_out_.shape()[1];
 
   // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
   double mu_ik  = mu_out_(i, k);
@@ -129,7 +129,7 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group,
 
 void
 IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
@@ -147,10 +147,10 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
   int n = energy_out_.shape()[1];
   if (!skewed_) {
     // All bins equally likely
-    j = prn(prn_seeds, stream) * n;
+    j = prn(prn_seed) * n;
   } else {
     // Distribution skewed away from edge points
-    double r = prn(prn_seeds, stream) * (n - 3);
+    double r = prn(prn_seed) * (n - 3);
     if (r > 1.0) {
       // equally likely N-4 middle bins
       j = r + 1;
@@ -178,7 +178,7 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
 
   // Sample outgoing cosine bin
   int m = mu_out_.shape()[2];
-  int k = prn(prn_seeds, stream) * m;
+  int k = prn(prn_seed) * m;
 
   // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
   double mu_ijk  = mu_out_(i, j, k);
@@ -233,7 +233,7 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
 
 void
 IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
@@ -241,7 +241,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
   get_energy_index(energy_, E_in, i, f);
 
   // Sample between ith and [i+1]th bin
-  int l = f > prn(prn_seeds, stream) ? i + 1 : i;
+  int l = f > prn(prn_seed) ? i + 1 : i;
 
   // Determine endpoints on grid i
   auto n = distribution_[i].e_out.size();
@@ -259,7 +259,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
   // Determine outgoing energy bin
   // (First reset n_energy_out to the right value)
   n = distribution_[l].n_e_out;
-  double r1 = prn(prn_seeds, stream);
+  double r1 = prn(prn_seed);
   double c_j = distribution_[l].e_out_cdf[0];
   double c_j1;
   std::size_t j;
@@ -298,7 +298,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
 
   // Sample outgoing cosine bin
   int n_mu = distribution_[l].mu.shape()[1];
-  std::size_t k = prn(prn_seeds, stream) * n_mu;
+  std::size_t k = prn(prn_seed) * 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
@@ -326,7 +326,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
   }
 
   // Smear angle
-  mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seeds, stream) - 0.5);
+  mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seed) - 0.5);
 }
 
 } // namespace openmc
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
index 1f6c1a9b7..6f65e0ef2 100644
--- a/src/secondary_uncorrelated.cpp
+++ b/src/secondary_uncorrelated.cpp
@@ -53,7 +53,7 @@ UncorrelatedAngleEnergy::UncorrelatedAngleEnergy(hid_t group)
 
 void
 UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
-  uint64_t* prn_seeds, int stream) const
+  uint64_t* prn_seed) const
 {
   // Sample cosine of scattering angle
   if (fission_) {
@@ -62,14 +62,14 @@ UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
     mu = 1.0;
     // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   } else if (!angle_.empty()) {
-    mu = angle_.sample(E_in, prn_seeds, stream);
+    mu = angle_.sample(E_in, prn_seed);
   } else {
     // no angle distribution given => assume isotropic for all energies
-    mu = 2.0*prn(prn_seeds, stream) - 1.0;
+    mu = 2.0*prn(prn_seed) - 1.0;
   }
 
   // Sample outgoing energy
-  E_out = energy_->sample(E_in, prn_seeds, stream);
+  E_out = energy_->sample(E_in, prn_seed);
 }
 
 } // namespace openmc
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 86f20a93d..405102827 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_);
+  init_particle_seeds(particle_seed, p->prn_seeds_);
 
   // set particle trace
   simulation::trace = false;
diff --git a/src/source.cpp b/src/source.cpp
index 320034735..483fd4b92 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -142,7 +142,7 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
 }
 
 
-Particle::Bank SourceDistribution::sample(uint64_t* prn_seeds, int stream) const
+Particle::Bank SourceDistribution::sample(uint64_t* prn_seed) const
 {
   Particle::Bank site;
 
@@ -158,7 +158,7 @@ Particle::Bank SourceDistribution::sample(uint64_t* prn_seeds, int stream) const
     site.particle = particle_;
 
     // Sample spatial distribution
-    site.r = space_->sample(prn_seeds, stream);
+    site.r = space_->sample(prn_seed);
     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(uint64_t* prn_seeds, int stream) const
   ++n_accept;
 
   // Sample angle
-  site.u = angle_->sample(prn_seeds, stream);
+  site.u = angle_->sample(prn_seed);
 
   // Check for monoenergetic source above maximum particle energy
   auto p = static_cast<int>(particle_);
@@ -218,7 +218,7 @@ Particle::Bank SourceDistribution::sample(uint64_t* prn_seeds, int stream) const
 
   while (true) {
     // Sample energy spectrum
-    site.E = energy_->sample(prn_seeds, stream);
+    site.E = energy_->sample(prn_seed);
 
     // 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,11 +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;
-      uint64_t prn_seeds[N_STREAMS];
-      set_particle_seed(id, prn_seeds);
+      uint64_t prn_seed;
+      init_seed(id, &prn_seed, STREAM_SOURCE);
 
       // sample external source distribution
-      simulation::source_bank[i] = sample_external_source(prn_seeds);
+      simulation::source_bank[i] = sample_external_source(&prn_seed);
     }
   }
 
@@ -288,11 +288,8 @@ void initialize_source()
   }
 }
 
-Particle::Bank sample_external_source(uint64_t* prn_seeds)
+Particle::Bank sample_external_source(uint64_t* prn_seed)
 {
-  // Set the random number generator to the source stream.
-  int stream = STREAM_SOURCE;
-
   // Determine total source strength
   double total_strength = 0.0;
   for (auto& s : model::external_sources)
@@ -301,7 +298,7 @@ Particle::Bank sample_external_source(uint64_t* prn_seeds)
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    double xi = prn(prn_seeds, stream)*total_strength;
+    double xi = prn(prn_seed)*total_strength;
     double c = 0.0;
     for (; i < model::external_sources.size(); ++i) {
       c += model::external_sources[i].strength();
@@ -310,7 +307,7 @@ Particle::Bank sample_external_source(uint64_t* prn_seeds)
   }
 
   // Sample source site from i-th source distribution
-  Particle::Bank site {model::external_sources[i].sample(prn_seeds, stream)};
+  Particle::Bank site {model::external_sources[i].sample(prn_seed)};
 
   // If running in MG, convert site.E to group
   if (!settings::run_CE) {
@@ -334,11 +331,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;
-      uint64_t prn_seeds[N_STREAMS];
-      set_particle_seed(id, prn_seeds);
+      uint64_t prn_seed;
+      init_seed(id, &prn_seed, STREAM_SOURCE);
 
       // sample external source distribution
-      simulation::source_bank[i] = sample_external_source(prn_seeds);
+      simulation::source_bank[i] = sample_external_source(&prn_seed);
     }
   }
 }
diff --git a/src/surface.cpp b/src/surface.cpp
index 958cd3851..feea34d8d 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, uint64_t* prn_seeds, int stream) const
+Surface::diffuse_reflect(Position r, Direction u, uint64_t* prn_seed) const
 {
   // Diffuse reflect direction according to the normal.
   // cosine distribution
@@ -208,10 +208,10 @@ Surface::diffuse_reflect(Position r, Direction u, uint64_t* prn_seeds, int strea
 
   // 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(prn_seeds, stream)) : std::sqrt(prn(prn_seeds, stream));
+                  -std::sqrt(prn(prn_seed)) : std::sqrt(prn(prn_seed));
 
   // sample azimuthal distribution uniformly
-  u = rotate_angle(n, mu, nullptr, prn_seeds, stream);
+  u = rotate_angle(n, mu, nullptr, prn_seed);
 
   // normalize the direction
   return u/u.norm();
diff --git a/src/thermal.cpp b/src/thermal.cpp
index ba0939e1e..93933e2c1 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -151,7 +151,7 @@ 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,
-                                uint64_t* prn_seeds, int stream) const
+                                uint64_t* prn_seed) const
 {
   // Determine temperature for S(a,b) table
   double kT = sqrtkT*sqrtkT;
@@ -173,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(prn_seeds, stream)) ++i;
+    if (f > prn(prn_seed)) ++i;
   }
 
   // Set temperature index
@@ -266,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, uint64_t* prn_seeds, int stream)
+                    double* E_out, double* mu, uint64_t* prn_seed)
 {
   // Determine whether inelastic or elastic scattering will occur
-  if (prn(prn_seeds, stream) < micro_xs.thermal_elastic / micro_xs.thermal) {
-    elastic_.distribution->sample(E, *E_out, *mu, prn_seeds, stream);
+  if (prn(prn_seed) < micro_xs.thermal_elastic / micro_xs.thermal) {
+    elastic_.distribution->sample(E, *E_out, *mu, prn_seed);
   } else {
-    inelastic_.distribution->sample(E, *E_out, *mu, prn_seeds, stream);
+    inelastic_.distribution->sample(E, *E_out, *mu, prn_seed);
   }
 
   // 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 0eb2df6f4..7df5047f6 100644
--- a/src/volume_calc.cpp
+++ b/src/volume_calc.cpp
@@ -125,16 +125,15 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
       std::vector<std::vector<int>> hits(n);
       Particle p;
 
-      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, prn_seeds);
+        uint64_t prn_seed;
+        int64_t id = iterations * n_samples_ + i;
+        init_seed(id, &prn_seed, STREAM_VOLUME);
 
         p.n_coord_ = 1;
-        Position xi {prn(prn_seeds, stream), prn(prn_seeds, stream), prn(prn_seeds, stream)};
+        Position xi {prn(&prn_seed), prn(&prn_seed), prn(&prn_seed)};
         p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
         p.u() = {0.5, 0.5, 0.5};
 
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 24dce4f95..455a1116e 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -156,13 +156,11 @@ def test_rotate_angle():
     uvw0 = np.array([1., 0., 0.])
     phi = 0.
     mu = 0.
-    prn_seeds = [1, 2, 3, 4, 5, 6]
-    stream = 0
 
     # reference: mu of 0 pulls the vector the bottom, so:
     ref_uvw = np.array([0., 0., -1.])
 
-    test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream)
+    test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi)
 
     assert np.array_equal(ref_uvw, test_uvw)
 
@@ -170,59 +168,56 @@ def test_rotate_angle():
     mu = 1.
     ref_uvw = np.array([1., 0., 0.])
 
-    test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream)
+    test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi)
 
     assert np.array_equal(ref_uvw, test_uvw)
 
     # Now to test phi is None
     mu = 0.9
     phi = None
+    prn_seed = 1
 
     # When seed = 1, phi will be sampled as 1.9116495709698769
     # The resultant reference is from hand-calculations given the above
     ref_uvw = [0.9, 0.410813051297112, 0.1457142302040]
-    test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream)
+    test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seed)
 
     assert np.allclose(ref_uvw, test_uvw)
 
 
 def test_maxwell_spectrum():
-    prn_seeds = [1, 2, 3, 4, 5, 6]
-    stream = 0
+    prn_seed = 1
     T = 0.5
     ref_val = 0.6129982175261098
-    test_val = openmc.lib.math.maxwell_spectrum(T, prn_seeds, stream)
+    test_val = openmc.lib.math.maxwell_spectrum(T, prn_seed)
 
     assert ref_val == test_val
 
 
 def test_watt_spectrum():
-    prn_seeds = [1, 2, 3, 4, 5, 6]
-    stream = 0
+    prn_seed = 1
     a = 0.5
     b = 0.75
     ref_val = 0.6247242713640233
-    test_val = openmc.lib.math.watt_spectrum(a, b, prn_seeds, stream)
+    test_val = openmc.lib.math.watt_spectrum(a, b, prn_seed)
 
     assert ref_val == test_val
 
 
 def test_normal_dist():
-    prn_seeds = [1, 2, 3, 4, 5, 6]
-    stream = 0
+    prn_seed = 1
     a = 14.08
     b = 0.0
     ref_val = 14.08
-    test_val = openmc.lib.math.normal_variate(a, b, prn_seeds, stream)
+    test_val = openmc.lib.math.normal_variate(a, b, prn_seed)
 
     assert ref_val == pytest.approx(test_val)
 
-    prn_seeds = [1, 2, 3, 4, 5, 6]
-    stream = 0
+    prn_seed = 1
     a = 14.08
     b = 1.0
     ref_val = 16.436645416691427
-    test_val = openmc.lib.math.normal_variate(a, b, prn_seeds, stream)
+    test_val = openmc.lib.math.normal_variate(a, b, prn_seed)
 
     assert ref_val == pytest.approx(test_val)