diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml
index 7b0a985f10..5845898c8d 100644
--- a/examples/xml/pincell/settings.xml
+++ b/examples/xml/pincell/settings.xml
@@ -3,8 +3,8 @@
 
   <!-- Define how many particles to run and for how many batches -->
   <run_mode>eigenvalue</run_mode>
-  <batches>10</batches>
-  <inactive>5</inactive>
+  <batches>100</batches>
+  <inactive>10</inactive>
   <particles>1000</particles>
 
   <!-- The starting source is a uniform distribution over the entire pin
diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h
index 2d1ebc8f3b..c88acd5696 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 c53511a739..d62ba1624a 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 c571aa61d0..23d7aa008e 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 8ff5b0fc04..dd59a23596 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 bc62386a08..08f1cd0b51 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 8d4aee5d5a..4895dd6ae3 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/geometry.h b/include/openmc/geometry.h
index a6e1e40631..e386348c44 100644
--- a/include/openmc/geometry.h
+++ b/include/openmc/geometry.h
@@ -28,14 +28,12 @@ extern std::vector<int64_t> overlap_check_count;
 // Information about nearest boundary crossing
 //==============================================================================
 
-/*
 struct BoundaryInfo {
   double distance {INFINITY};   //!< distance to nearest boundary
   int surface_index {0}; //!< if boundary is surface, index in surfaces vector
   int coord_level;   //!< coordinate level after crossing boundary
   std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
 };
-*/
 
 //==============================================================================
 //! Check two distances by coincidence tolerance
diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index e2642c1a33..1440d13b9b 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 st
 //!
 //! @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_see
 //! @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 cfb696a62e..131d248e11 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/particle.h b/include/openmc/particle.h
index 3a9f8ff163..e99aaf9fc0 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -129,21 +129,8 @@ struct MacroXS {
   double incoherent;      //!< macroscopic incoherent xs
   double photoelectric;   //!< macroscopic photoelectric xs
   double pair_production; //!< macroscopic pair production xs
-  int i_grid;
 };
 
-//==============================================================================
-// Information about nearest boundary crossing
-//==============================================================================
-
-struct BoundaryInfo {
-  double distance {INFINITY};   //!< distance to nearest boundary
-  int surface_index {0}; //!< if boundary is surface, index in surfaces vector
-  int coord_level;   //!< coordinate level after crossing boundary
-  std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
-};
-
-
 //============================================================================
 //! State of a particle being transported through geometry
 //============================================================================
@@ -288,7 +275,6 @@ public:
 
   // Indices for various arrays
   int surface_ {0};             //!< index for surface particle is on
-  BoundaryInfo boundary_;
   int cell_born_ {-1};      //!< index for cell particle was born in
   int material_ {-1};       //!< index for current material
   int material_last_ {-1};  //!< index for last material
@@ -313,14 +299,4 @@ public:
 
 } // namespace openmc
 
-extern template class std::vector<openmc::Particle::Bank>;
-
-namespace openmc{
-
-extern std::vector<Particle> particle_bank;
-#pragma omp threadprivate(particle_bank)
-
-}
-
-
 #endif // OPENMC_PARTICLE_H
diff --git a/include/openmc/photon.h b/include/openmc/photon.h
index f5b79455ce..c7bbeddd2e 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 6726e56dc4..fce164cd6e 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 b99f914a4e..350fa7d245 100644
--- a/include/openmc/random_lcg.h
+++ b/include/openmc/random_lcg.h
@@ -10,53 +10,64 @@ namespace openmc {
 // Module constants.
 //==============================================================================
 
-constexpr int N_STREAMS = 6;
-extern "C" const int STREAM_TRACKING;
-extern "C" const int STREAM_TALLIES;
-extern "C" const int STREAM_SOURCE;
-extern "C" const int STREAM_URR_PTABLE;
-extern "C" const int STREAM_VOLUME;
-extern "C" const int STREAM_PHOTON;
-constexpr int64_t DEFAULT_SEED = 1;
+constexpr int N_STREAMS         {6};
+constexpr int STREAM_TRACKING   {0};
+constexpr int STREAM_TALLIES    {1};
+constexpr int STREAM_SOURCE     {2};
+constexpr int STREAM_URR_PTABLE {3};
+constexpr int STREAM_VOLUME     {4};
+constexpr int STREAM_PHOTON     {5};
+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);
+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);
+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 );
+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
+//==============================================================================
+
+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);
+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 23a4799f3d..c61eec6b01 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 458a4c24d7..8f3a79b77d 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 595fc43911..6b7c367961 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 90ad1e7102..b9f1906353 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 d8a1dafdf5..b810d17d01 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 edd7356f62..74bb3265f2 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 5cde70b4a3..66ef834f19 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/simulation.h b/include/openmc/simulation.h
index 7981286be2..8efac4c7f8 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -16,10 +16,6 @@ constexpr int STATUS_EXIT_NORMAL {0};
 constexpr int STATUS_EXIT_MAX_BATCH {1};
 constexpr int STATUS_EXIT_ON_TRIGGER {2};
 
-extern Particle::Bank * shared_fission_bank;
-extern int shared_fission_bank_length;
-extern int shared_fission_bank_max;
-
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -42,9 +38,7 @@ extern "C" int restart_batch;   //!< batch at which a restart job resumed
 extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
 extern "C" int total_gen;        //!< total number of generations simulated
 extern double total_weight;  //!< Total source weight in a batch
-extern int64_t thread_work_index;
-extern int64_t work_per_rank;       //!< number of particles per MPI rank
-extern int64_t work_per_thread;    //!< number of particles on a given tread
+extern int64_t work_per_rank;         //!< number of particles per MPI rank
 
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;
@@ -55,8 +49,7 @@ extern std::vector<int64_t> work_index;
 // Threadprivate variables
 extern "C" bool trace;     //!< flag to show debug information
 
-
-#pragma omp threadprivate(current_work, work_per_thread, trace)
+#pragma omp threadprivate(current_work, trace)
 
 } // namespace simulation
 
diff --git a/include/openmc/source.h b/include/openmc/source.h
index 5826fe93e7..8a7c30a09d 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 8060afcd1f..35fa84bd13 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 b0a50de657..37ad8f4cbd 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 1d0265e039..991df8f9cc 100644
--- a/openmc/lib/math.py
+++ b/openmc/lib/math.py
@@ -1,10 +1,12 @@
-from ctypes import (c_int, c_double, POINTER, c_uint64)
+from ctypes import c_int, c_double, POINTER, c_uint64
 
 import numpy as np
 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 a99d118f0c..85733ecccb 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 0b56f9857e..e48ce8baae 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,15 +108,15 @@ 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);
 }
 
 //==============================================================================
 // Normal implementation
 //==============================================================================
-Normal::Normal(pugi::xml_node node) 
+Normal::Normal(pugi::xml_node node)
 {
   auto params = get_node_array<double>(node,"parameters");
   if (params.size() != 2)
@@ -127,15 +127,15 @@ 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);
 }
 
 //==============================================================================
 // Muir implementation
 //==============================================================================
-Muir::Muir(pugi::xml_node node) 
+Muir::Muir(pugi::xml_node node)
 {
   auto params = get_node_array<double>(node,"parameters");
   if (params.size() != 3)
@@ -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 f946d0602e..076f824701 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,11 +83,10 @@ double AngleDistribution::sample(double E, uint64_t * prn_seeds, int stream) con
   }
 
   // Sample between the ith and (i+1)th bin
-  if (r > prn(prn_seeds, stream)) ++i;
+  if (r > prn(prn_seed)) ++i;
 
-  //assert(i >= 0 );
   // 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 48b86b8cc7..3e19951451 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) con
   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) con
   // 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 df2f8d682d..7112342018 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 9dd264107a..ce2b7f34a2 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 747003d9ba..cd0699a32b 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -114,17 +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.
 
-  //std::cout << "MPI rank " << mpi::rank << " start = " << start << std::endl;
-  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.
@@ -159,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 6e639432b3..21afd696c3 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 9dceb5127c..e2195d512c 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, i
   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, i
 //==============================================================================
 
 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 b5da3f60aa..bc0fed2e0f 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/output.cpp b/src/output.cpp
index 29d60830eb..899bce86e2 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -67,15 +67,7 @@ void title()
     "                     ###############     %%%%%%%%%%%%%%%%\n" <<
     "                       ############     %%%%%%%%%%%%%%%\n" <<
     "                          ########     %%%%%%%%%%%%%%\n" <<
-    "                                      %%%%%%%%%%%\n" <<
-    "                              )                             (        (\n" << 
-    "                           ( /(   *   )        (     (      )\\ )     )\\ )\n" <<  
-    "         (    (   (   (    )\\())` )  /(      ( )\\    )\\    (()/( (  (()/(\n" << 
-    "         )\\   )\\  )\\  )\\  ((_)\\  ( )(_))___  )((_)((((_)(   /(_)))\\  /(_))\n" << 
-    "        ((_) ((_)((_)((_)  _((_)(_(_())|___|((_)_  )\\ _ )\\ (_)) ((_)(_))_\n" <<  
-    "        | __|\\ \\ / / | __|| \\| ||_   _|      | _ ) (_)_\\(_)/ __|| __||   \\\n" << 
-    "        | _|  \\ V /  | _| | .` |  | |        | _ \\  / _ \\  \\__ \\| _| | |) |\n" << 
-    "        |___|  \\_/   |___||_|\\_|  |_|        |___/ /_/ \\_\\ |___/|___||___/\n\n";
+    "                                      %%%%%%%%%%%\n\n";
 
   // Write version information
   std::cout <<
diff --git a/src/particle.cpp b/src/particle.cpp
index 662c2b342f..eb5d43b27e 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -27,14 +27,8 @@
 #include "openmc/tallies/tally_scoring.h"
 #include "openmc/track_output.h"
 
-// Explicit template instantiation definition
-template class std::vector<openmc::Particle>;
-
-
 namespace openmc {
 
-std::vector<Particle> particle_bank;
-
 //==============================================================================
 // LocalCoord implementation
 //==============================================================================
@@ -140,15 +134,14 @@ void
 Particle::transport()
 {
   // Display message if high verbosity or trace is on
-  //if (settings::verbosity >= 9 || simulation::trace) {
-    // write_message("Simulating Particle " + std::to_string(id_));
-  //}
+  if (settings::verbosity >= 9 || simulation::trace) {
+     write_message("Simulating Particle " + std::to_string(id_));
+  }
 
   // Initialize number of events to zero
   int n_event = 0;
 
   // Add paricle's starting weight to count for normalizing tallies later
-  /*
   #pragma omp atomic
   simulation::total_weight += wgt_;
 
@@ -162,7 +155,6 @@ Particle::transport()
 
   // Every particle starts with no accumulated flux derivative.
   if (!model::active_tallies.empty()) zero_flux_derivs();
-  */
 
   while (true) {
     // Set the random number stream
@@ -238,7 +230,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
@@ -474,7 +466,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 17eb5a3e6b..b7d94dd807 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 803d891846..6cf43cbc52 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, i
     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, i
     // 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
     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
   } 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 893519993d..7fad84af67 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -79,7 +79,6 @@ void collision(Particle* p)
 
 void sample_neutron_reaction(Particle* p)
 {
-	//std::cout << "particle energy e = " << p->E_ << std::endl;
   // Sample a nuclide within the material
   int i_nuclide = sample_nuclide(p);
 
@@ -130,18 +129,15 @@ void sample_neutron_reaction(Particle* p)
   }
   if (!p->alive_) return;
 
-	//std::cout << "before scatter particle energy e = " << p->E_ << " with nuclide = " << i_nuclide <<  std::endl;
   // Sample a scattering reaction and determine the secondary energy of the
   // exiting neutron
   scatter(p, i_nuclide);
-	
-  //std::cout << "after scatter particle energy e = " << p->E_ << std::endl;
 
   // Advance URR seed stream 'N' times after energy changes
   if (p->E_ != p->E_last_) {
-	p->stream_ = STREAM_URR_PTABLE;
-    advance_prn_seed(data::nuclides.size(), p->prn_seeds_, p->stream_);
-	p->stream_ = STREAM_TRACKING;
+    p->stream_ = STREAM_URR_PTABLE;
+    advance_prn_seed(data::nuclides.size(), p->prn_seeds_ + p->stream_);
+    p->stream_ = STREAM_TRACKING;
   }
 
   // Play russian roulette if survival biasing is turned on
@@ -165,7 +161,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
@@ -177,16 +173,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 
   p->fission_ = true;
   for (int i = 0; i < nu; ++i) {
-	  /*
     // Create new bank site and get reference to last element
-    bank.emplace_back();
-    auto& site {bank.back()};
-
-    // Bank source neutrons by copying the particle data
-    site.r = p->r();
-    site.particle = Particle::Type::neutron;
-    site.wgt = 1. / weight;
-	*/
 	  int idx;
 	  #pragma omp atomic capture
 	  idx = shared_fission_bank_length++;
@@ -197,12 +184,10 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 	  site->parent_id = p->id_;
 
     // 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;
-    p->delayed_group_ = site->delayed_group;
+    p->delayed_group_ = site.delayed_group;
 
     // Increment the number of neutrons born delayed
     if (p->delayed_group_ > 0) {
@@ -241,13 +226,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;
@@ -258,7 +243,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.
@@ -270,13 +255,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);
     }
 
@@ -290,7 +275,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;
@@ -322,8 +307,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);
@@ -331,7 +316,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);
@@ -359,14 +344,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;
@@ -400,8 +385,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);
@@ -419,7 +404,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_]};
@@ -444,7 +429,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_]};
@@ -497,7 +482,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
@@ -524,7 +509,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
@@ -572,7 +557,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_[
@@ -600,7 +585,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
@@ -616,7 +601,6 @@ void scatter(Particle* p, int i_nuclide)
     // Determine temperature
     double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp];
 
-	//std::cout << "we are doing an elastic scatter" << std::endl;
     // Perform collision physics for elastic scattering
     elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p);
 
@@ -629,7 +613,6 @@ void scatter(Particle* p, int i_nuclide)
     // =======================================================================
     // S(A,B) SCATTERING
 
-	//std::cout << "we are doing an SAB scatter" << std::endl;
     sab_scatter(i_nuclide, micro.index_sab, p);
 
     p->event_mt_ = ELASTIC;
@@ -637,7 +620,6 @@ void scatter(Particle* p, int i_nuclide)
   }
 
   if (!sampled) {
-	//std::cout << "we are doing an Inelastic scatter" << std::endl;
     // =======================================================================
     // INELASTIC SCATTERING
 
@@ -664,9 +646,7 @@ void scatter(Particle* p, int i_nuclide)
 
     // Perform collision physics for inelastic scattering
     const auto& rx {nuc->reactions_[i]};
-	//std::cout << "Energy before Inelastic scatter = " << p->E_ << std::endl;
     inelastic_scatter(nuc.get(), rx.get(), p);
-	//std::cout << "Energy after Inelastic scatter = " << p->E_ << std::endl;
     p->event_mt_ = rx->mt_;
   }
 
@@ -679,8 +659,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;
@@ -707,7 +687,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
@@ -725,9 +705,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
@@ -736,7 +716,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;
@@ -765,15 +745,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;
@@ -805,7 +785,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: {
@@ -839,7 +819,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) {
@@ -863,7 +843,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;
@@ -872,7 +852,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) {
@@ -892,10 +872,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];
@@ -913,7 +893,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);
         }
       }
     }
@@ -924,7 +904,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);
@@ -933,10 +913,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
@@ -948,7 +928,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;
     }
 
@@ -956,14 +936,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
@@ -971,19 +951,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);
@@ -994,12 +974,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) {
@@ -1023,7 +1003,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);
@@ -1048,7 +1028,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);
@@ -1073,40 +1053,33 @@ 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
   if (rx->scatter_in_cm_) {
     double E_cm = E;
-	//std::cout << "E_cm = " << E_cm << std::endl;
 
     // determine outgoing energy in lab
     double A = nuc->awr_;
-//	std::cout << "A = " << nuc->awr_ << std::endl;
     E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm))
           / ((A + 1.0)*(A + 1.0));
-	// here's where it goes wrong
-	//std::cout << "E = " << E << std::endl;
 
     // determine outgoing angle in lab
     mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E);
-//	std::cout << "mu = " << mu << std::endl;
   }
 
   // Because of floating-point roundoff, it may be possible for mu to be
   // outside of the range [-1,1). In these cases, we just set mu to exactly -1
   // or 1
   if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
-//	std::cout << "mu = " << mu << std::endl;
 
   // Set outgoing energy and scattering angle
-//	std::cout << "FINAL E = " << E << " mu = " << mu << std::endl;
   p->E_ = E;
   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);
@@ -1127,7 +1100,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) {
@@ -1140,10 +1113,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 e5b43e2234..c26ebdfd2c 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 00e2fd2611..6f198f1e56 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -79,10 +79,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_];
@@ -104,7 +104,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++;
   }
 
@@ -129,10 +129,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);
@@ -140,8 +140,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
@@ -181,7 +181,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 07c6129520..dbc1ccb9e0 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 6c76cc6cfb..851d72e3e2 100644
--- a/src/random_lcg.cpp
+++ b/src/random_lcg.cpp
@@ -6,15 +6,6 @@
 
 namespace openmc {
 
-
-// Constants
-extern "C" const int STREAM_TRACKING   {0};
-extern "C" const int STREAM_TALLIES    {1};
-extern "C" const int STREAM_SOURCE     {2};
-extern "C" const int STREAM_URR_PTABLE {3};
-extern "C" const int STREAM_VOLUME     {4};
-extern "C" const int STREAM_PHOTON     {5};
-
 // Starting seed
 int64_t master_seed {1};
 
@@ -32,34 +23,40 @@ constexpr double   prn_norm   {1.0 / prn_mod};           // 2^-63
 // PRN
 //==============================================================================
 
-extern "C" double
-prn(uint64_t * prn_seeds, int stream)
+double 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;
 }
 
 //==============================================================================
 // FUTURE_PRN
 //==============================================================================
 
-extern "C" double
-future_prn(int64_t n, uint64_t * prn_seeds, int stream)
+double 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
+// INIT_SEED
 //==============================================================================
 
-extern "C" void
-set_particle_seed(int64_t id, uint64_t * prn_seeds)
+void 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);
+}
+
+//==============================================================================
+// INIT_PARTICLE_SEEDS
+//==============================================================================
+
+void init_particle_seeds(int64_t id, uint64_t* prn_seeds)
 {
 	//std::cout << "Master seed = " << master_seed << " id = " << id << std::endl;
   for (int i = 0; i < N_STREAMS; i++) {
@@ -71,18 +68,16 @@ set_particle_seed(int64_t id, uint64_t * prn_seeds)
 // ADVANCE_PRN_SEED
 //==============================================================================
 
-extern "C" void
-advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream)
+void 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);
 }
 
 //==============================================================================
 // FUTURE_SEED
 //==============================================================================
 
-uint64_t
-future_seed(uint64_t n, uint64_t seed)
+uint64_t future_seed(uint64_t n, uint64_t prn_seed)
 {
   // Make sure nskip is less than 2^M.
   n &= prn_mask;
@@ -113,7 +108,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;
 }
 
 //==============================================================================
@@ -122,8 +117,7 @@ future_seed(uint64_t n, uint64_t seed)
 
 extern "C" int64_t openmc_get_seed() {return master_seed;}
 
-extern "C" void
-openmc_set_seed(int64_t new_seed)
+extern "C" void openmc_set_seed(int64_t new_seed)
 {
   master_seed = new_seed;
 }
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
index e047be8de4..d92d14517a 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 a5eeae6211..5d06375e50 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 4c3b93262e..81e1ec1f9f 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 781923720d..b0f9854ac7 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_
   }
 
   // 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_
   // 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_
   }
 
   // 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 5d32d3e860..46e56145b2 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 fc4df50a6a..13b5746fc6 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 f4bddeebeb..6f65e0ef26 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 6b9e91ed81..ded61b9833 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -2,32 +2,23 @@
 
 #include "openmc/bank.h"
 #include "openmc/capi.h"
-#include "openmc/cell.h"
 #include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
-#include "openmc/geometry.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
-#include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
 #include "openmc/particle.h"
 #include "openmc/photon.h"
-#include "openmc/physics.h"
- #include "openmc/physics_mg.h"
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
 #include "openmc/source.h"
 #include "openmc/state_point.h"
-#include "openmc/thermal.h"
 #include "openmc/timer.h"
-#include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
-#include "openmc/tallies/tally_scoring.h"
 #include "openmc/tallies/trigger.h"
-#include "openmc/track_output.h"
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -42,6 +33,7 @@
 #include <string>
 #include <chrono> 
 
+
 namespace openmc {
 
 	/*
@@ -791,12 +783,10 @@ int openmc_run()
 
   int err = 0;
   int status = 0;
-
   while (status == 0 && err == 0) {
     err = openmc_next_batch(&status);
   }
 
-
   openmc_simulation_finalize();
   return err;
 }
@@ -807,16 +797,9 @@ int openmc_simulation_init()
 
   // Skip if simulation has already been initialized
   if (simulation::initialized) return 0;
-	
 
   // Determine how much work each process should do
   calculate_work();
-  
-  // We assume that the number of fission sites will be at most
-	// 10x the number of neutrons run on this rank
-	int shared_fission_bank_size = 10 * simulation::work_per_rank;
-	// Initializes the shared fission bank
-	init_shared_fission_bank(shared_fission_bank_size);
 
   // Allocate array for matching filter bins
   #pragma omp parallel
@@ -829,11 +812,10 @@ int openmc_simulation_init()
   allocate_banks();
 
   // Allocate tally results arrays if they're not allocated yet
-	  //std::cout << "Trying to Initialize Results..." << std::endl;
+
   for (auto& t : model::tallies) {
     t->init_results();
   }
-	  //std::cout << "Success in Initializing Results!" << std::endl;
 
   // Set up material nuclide index mapping
   for (auto& mat : model::materials) {
@@ -878,8 +860,6 @@ int openmc_simulation_finalize()
 
   // Skip if simulation was never run
   if (!simulation::initialized) return 0;
-  
-  free_shared_fission_bank();
 
   // Stop active batch timer and start finalization timer
   simulation::time_active.stop();
@@ -947,29 +927,9 @@ int openmc_next_batch(int* status)
     // Start timer for transport
     simulation::time_transport.start();
 
-    // ====================================================================
-    // LOOP OVER PARTICLES
+	  simulation::current_work = 1;
 
-	simulation::current_work = 1;
-
-    //#pragma omp parallel
-    {
-      transport();
-    }
-
-	/*
-    #pragma omp parallel for schedule(runtime)
-    for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
-      simulation::current_work = i_work;
-
-      // grab source particle from bank
-      Particle p;
-      initialize_history(&p, simulation::current_work);
-
-      // transport particle
-      p.transport();
-    }
-	*/
+    transport();
 
     // Accumulate time for transport
     simulation::time_transport.stop();
@@ -1026,9 +986,7 @@ int restart_batch;
 bool satisfy_triggers {false};
 int total_gen {0};
 double total_weight;
-int64_t thread_work_index;
 int64_t work_per_rank;
-int64_t work_per_thread;
 
 const RegularMesh* entropy_mesh {nullptr};
 const RegularMesh* ufs_mesh {nullptr};
@@ -1252,30 +1210,14 @@ void finalize_generation()
   }
 
 
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-//#ifdef _OPENMP
-    // Join the fission bank from each thread into one global fission bank
-    //join_bank_from_threads();
-	
+  if (settings::run_mode == RUN_MODE_EIGENVALUE) {	
 	  // We need to move all the stuff from the shared_fission_bank into the real one.
-		//std::vector<Particle::Bank> shared_fission_bank_vector(shared_fission_bank, shared_fission_bank + shared_fission_bank_length);
-      //simulation::fission_bank = shared_fission_bank_vector;
-	  //std::cout << "Fission bank length = " << shared_fission_bank_length << std::endl;
 	  for( int i = 0; i < shared_fission_bank_length; i++ )
 		  simulation::fission_bank.push_back(shared_fission_bank[i]);
 	  shared_fission_bank_length = 0;
-//#endif
-      // Sorts the fission bank so as to allow for reproducibility
-      //std::sort(simulation::fission_bank.begin(), simulation::fission_bank.end(), bank_site_comparator());
-      //std::sort(simulation::fission_bank.begin(), simulation::fission_bank.end());
-      std::stable_sort(simulation::fission_bank.begin(), simulation::fission_bank.end());
-	  /*
-	  std::cout << "Fission bank on rank " << mpi::rank << " is of length " << simulation::fission_bank.size() << std::endl;
-	  for (Particle::Bank p : simulation::fission_bank )
-	  {
-	   std::cout <<  "E = " << p.E << std::endl;
-	  }
-	  */
+    
+    // Sorts the fission bank so as to allow for reproducibility
+    std::stable_sort(simulation::fission_bank.begin(), simulation::fission_bank.end());
 
     // Distribute fission bank across processors evenly
     synchronize_bank();
@@ -1302,10 +1244,8 @@ void finalize_generation()
 
 void initialize_history(Particle* p, int64_t index_source)
 {
-	//assert(index_source - 1 >= 0 );
   // set defaults
   p->from_source(&simulation::source_bank[index_source - 1]);
-  //p->from_source(&simulation::source_bank[index_source]);
 
   // set identifier for particle
   p->id_ = simulation::work_index[mpi::rank] + index_source;
@@ -1313,8 +1253,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_;
-  //std::cout << "MPI rank " << mpi::rank << " is initializing particle " << index_source << " with ID/seed " << particle_seed << std::endl;
-  set_particle_seed(particle_seed, p->prn_seeds_);
+  init_particle_seeds(particle_seed, p->prn_seeds_);
 
   // set particle trace
   simulation::trace = false;
@@ -1337,6 +1276,7 @@ void initialize_history(Particle* p, int64_t index_source)
     }
   }
 
+
 // Display message if high verbosity or trace is on
    if (settings::verbosity >= 9 || simulation::trace) {
       write_message("Simulating Particle " + std::to_string(p->id_));
@@ -1396,32 +1336,6 @@ void calculate_work()
     i_bank += work_i;
     simulation::work_index[i + 1] = i_bank;
   }
-	#ifdef _OPENMP
-   // Determine work per thread
-   int remainder_thread = simulation::work_per_rank % omp_get_max_threads();
-
-    #pragma omp parallel
-   {
-     simulation::work_per_thread = simulation::work_per_rank / omp_get_num_threads();
-     if (omp_get_thread_num() < remainder_thread) {
-       ++simulation::work_per_thread;
-     }
-   }
-
-    int64_t work_i = 0;
-   #pragma omp parallel for ordered
-   for (int i = 0; i < omp_get_num_threads(); ++i) {
-     #pragma omp ordered
-     {
-       simulation::thread_work_index = work_i;
-       work_i += simulation::work_per_thread;
-	   /*
-       std::cout << "Thread " << omp_get_thread_num() << ": " << simulation::work_per_thread
-         << std::endl;
-		 */
-     }
-   }
- #endif
 }
 
 #ifdef OPENMC_MPI
diff --git a/src/source.cpp b/src/source.cpp
index 860584a2c2..483fd4b92d 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) cons
     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) cons
   ++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) cons
 
   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 7c18d1c342..feea34d8d6 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -197,24 +197,24 @@ 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 
-  
+  // cosine distribution
+
   Direction n = this->normal(r);
   n /= n.norm();
   const double projection = n.dot(u);
-  
-  // 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));  
-  
-  // sample azimuthal distribution uniformly 
-  u = rotate_angle(n, mu, nullptr, prn_seeds, stream);
-  
-  // normalize the direction 
-  return u/u.norm();    
+
+  // 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_seed)) : std::sqrt(prn(prn_seed));
+
+  // sample azimuthal distribution uniformly
+  u = rotate_angle(n, mu, nullptr, prn_seed);
+
+  // normalize the direction
+  return u/u.norm();
 }
 
 CSGSurface::CSGSurface() : Surface{} {};
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 21474dc652..852918a0cd 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -804,7 +804,6 @@ Tally::init_triggers(pugi::xml_node node)
 void Tally::init_results()
 {
   int n_scores = scores_.size() * nuclides_.size();
-  //std::cout << "TRYING TO ALLOCATE n filter bins, n scores =  " << n_filter_bins_ << " " <<  n_scores << std::endl;
   results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
 }
 
diff --git a/src/thermal.cpp b/src/thermal.cpp
index 0af5b1614a..93933e2c10 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 0eb2df6f4c..7df5047f6b 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 24dce4f950..455a1116e4 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)