Merge pull request #1419 from jtramm/particle_owned_seeds

Particle owned random number state variables

include/openmc/angle_energy.h

@@ -1,6 +1,8 @@
 #ifndef OPENMC_ANGLE_ENERGY_H
 #define OPENMC_ANGLE_ENERGY_H
 
+#include <cstdint>
+
 namespace openmc {
 
 //==============================================================================
@@ -12,7 +14,8 @@ namespace openmc {
 
 class AngleEnergy {
 public:
-  virtual void sample(double E_in, double& E_out, double& mu) const = 0;
+  virtual void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const = 0;
   virtual ~AngleEnergy() = default;
 };
 

include/openmc/distribution.h

@@ -21,7 +21,7 @@ namespace openmc {
 class Distribution {
 public:
   virtual ~Distribution() = default;
-  virtual double sample() const = 0;
+  virtual double sample(uint64_t* seed) const = 0;
 };
 
 //==============================================================================
@@ -34,8 +34,9 @@ public:
   Discrete(const double* x, const double* p, int n);
 
   //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 
   // Properties
   const std::vector<double>& x() const { return x_; }
@@ -58,8 +59,9 @@ public:
   Uniform(double a, double b) : a_{a}, b_{b} {};
 
   //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 private:
   double a_; //!< Lower bound of distribution
   double b_; //!< Upper bound of distribution
@@ -75,8 +77,9 @@ public:
   Maxwell(double theta) : theta_{theta} { };
 
   //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 private:
   double theta_; //!< Factor in exponential [eV]
 };
@@ -91,8 +94,9 @@ public:
   Watt(double a, double b) : a_{a}, b_{b} { };
 
   //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 private:
   double a_; //!< Factor in exponential [eV]
   double b_; //!< Factor in square root [1/eV]
@@ -108,8 +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 seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 private:
   double mean_value_;    //!< middle of distribution [eV]
   double std_dev_; //!< standard deviation [eV]
@@ -126,8 +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 seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 private:
   // example DT fusion m_rat = 5 (D = 2 + T = 3)
   // ion temp = 20000 eV
@@ -148,8 +154,9 @@ public:
           const double* c=nullptr);
 
   //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 
   // x property
   std::vector<double>& x() { return x_; }
@@ -178,8 +185,9 @@ public:
   Equiprobable(const double* x, int n) : x_{x, x+n} { };
 
   //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
-  double sample() const;
+  double sample(uint64_t* seed) const;
 private:
   std::vector<double> x_; //! Possible outcomes
 };

include/openmc/distribution_angle.h

@@ -23,8 +23,9 @@ public:
 
   //! Sample an angle given an incident particle energy
   //! \param[in] E Particle energy in [eV]
+  //! \param[inout] seed pseudorandom number seed pointer
   //! \return Cosine of the angle in the range [-1,1]
-  double sample(double E) const;
+  double sample(double E, uint64_t* seed) const;
 
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty

include/openmc/distribution_energy.h

@@ -22,7 +22,7 @@ namespace openmc {
 
 class EnergyDistribution {
 public:
-  virtual double sample(double E) const = 0;
+  virtual double sample(double E, uint64_t* seed) const = 0;
   virtual ~EnergyDistribution() = default;
 };
 
@@ -36,8 +36,9 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E) const;
+  double sample(double E, uint64_t* seed) const;
 private:
   int primary_flag_; //!< Indicator of whether the photon is a primary or
                      //!< non-primary photon.
@@ -55,8 +56,9 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E) const;
+  double sample(double E, uint64_t* seed) const;
 private:
   double threshold_; //!< Energy threshold in lab, (A + 1)/A * |Q|
   double mass_ratio_; //!< (A/(A+1))^2
@@ -74,8 +76,9 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E) const;
+  double sample(double E, uint64_t* seed) const;
 private:
   //! Outgoing energy for a single incoming energy
   struct CTTable {
@@ -103,8 +106,9 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E) const;
+  double sample(double E, uint64_t* seed) const;
 private:
   Tabulated1D theta_; //!< Incoming energy dependent parameter
   double u_; //!< Restriction energy
@@ -121,8 +125,9 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E) const;
+  double sample(double E, uint64_t* seed) const;
 private:
   Tabulated1D theta_; //!< Incoming energy dependent parameter
   double u_; //!< Restriction energy
@@ -139,8 +144,9 @@ public:
 
   //! Sample energy distribution
   //! \param[in] E Incident particle energy in [eV]
+  //! \param[inout] seed Pseudorandom number seed pointer
   //! \return Sampled energy in [eV]
-  double sample(double E) const;
+  double sample(double E, uint64_t* seed) const;
 private:
   Tabulated1D a_; //!< Energy-dependent 'a' parameter
   Tabulated1D b_; //!< Energy-dependent 'b' parameter

include/openmc/distribution_multi.h

@@ -23,8 +23,9 @@ public:
   virtual ~UnitSphereDistribution() = default;
 
   //! Sample a direction from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Direction sampled
-  virtual Direction sample() const = 0;
+  virtual Direction sample(uint64_t* seed) const = 0;
 
   Direction u_ref_ {0.0, 0.0, 1.0};  //!< reference direction
 };
@@ -39,8 +40,9 @@ public:
   explicit PolarAzimuthal(pugi::xml_node node);
 
   //! Sample a direction from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Direction sampled
-  Direction sample() const;
+  Direction sample(uint64_t* seed) const;
 private:
   UPtrDist mu_;  //!< Distribution of polar angle
   UPtrDist phi_; //!< Distribution of azimuthal angle
@@ -55,8 +57,9 @@ public:
   Isotropic() { };
 
   //! Sample a direction from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled direction
-  Direction sample() const;
+  Direction sample(uint64_t* seed) const;
 };
 
 //==============================================================================
@@ -69,8 +72,9 @@ public:
   explicit Monodirectional(pugi::xml_node node) : UnitSphereDistribution{node} { };
 
   //! Sample a direction from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled direction
-  Direction sample() const;
+  Direction sample(uint64_t* seed) const;
 };
 
 using UPtrAngle = std::unique_ptr<UnitSphereDistribution>;

include/openmc/distribution_spatial.h

@@ -17,7 +17,7 @@ public:
   virtual ~SpatialDistribution() = default;
 
   //! Sample a position from the distribution
-  virtual Position sample() const = 0;
+  virtual Position sample(uint64_t* seed) const = 0;
 };
 
 //==============================================================================
@@ -29,8 +29,9 @@ public:
   explicit CartesianIndependent(pugi::xml_node node);
 
   //! Sample a position from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample() const;
+  Position sample(uint64_t* seed) const;
 private:
   UPtrDist x_; //!< Distribution of x coordinates
   UPtrDist y_; //!< Distribution of y coordinates
@@ -46,8 +47,9 @@ public:
   explicit SphericalIndependent(pugi::xml_node node);
 
   //! Sample a position from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample() const;
+  Position sample(uint64_t* seed) const;
 private:
   UPtrDist r_; //!< Distribution of r coordinates
   UPtrDist theta_; //!< Distribution of theta coordinates
@@ -64,8 +66,9 @@ public:
   explicit SpatialBox(pugi::xml_node node, bool fission=false);
 
   //! Sample a position from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample() const;
+  Position sample(uint64_t* seed) const;
 
   // Properties
   bool only_fissionable() const { return only_fissionable_; }
@@ -86,8 +89,9 @@ public:
   explicit SpatialPoint(pugi::xml_node node);
 
   //! Sample a position from the distribution
+  //! \param seed Pseudorandom number seed pointer
   //! \return Sampled position
-  Position sample() const;
+  Position sample(uint64_t* seed) const;
 private:
   Position r_; //!< Single position at which sites are generated
 };

include/openmc/math_functions.h

@@ -129,11 +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 seed A pointer to the pseudorandom seed
 //==============================================================================
 
-extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi);
+extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi,
+  uint64_t* seed);
 
-Direction rotate_angle(Direction u, double mu, const double* phi);
+Direction rotate_angle(Direction u, double mu, const double* phi,
+  uint64_t* seed);
 
 //==============================================================================
 //! Samples an energy from the Maxwell fission distribution based on a direct
@@ -144,10 +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 seed A pointer to the pseudorandom seed
 //! \return The sampled outgoing energy
 //==============================================================================
 
-extern "C" double maxwell_spectrum(double T);
+extern "C" double maxwell_spectrum(double T, uint64_t* seed);
 
 //==============================================================================
 //! Samples an energy from a Watt energy-dependent fission distribution.
@@ -159,10 +163,11 @@ extern "C" double maxwell_spectrum(double T);
 //!
 //! \param a Watt parameter a
 //! \param b Watt parameter b
+//! \param seed A pointer to the pseudorandom seed
 //! \return The sampled outgoing energy
 //==============================================================================
 
-extern "C" double watt_spectrum(double a, double b);
+extern "C" double watt_spectrum(double a, double b, uint64_t* seed);
 
 //==============================================================================
 //! Samples an energy from the Gaussian energy-dependent fission distribution.
@@ -175,10 +180,11 @@ extern "C" double watt_spectrum(double a, double b);
 //!
 //! @param mean mean of the Gaussian distribution
 //! @param std_dev standard deviation of the Gaussian distribution
+//! @param seed A pointer to the pseudorandom seed
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double normal_variate(double mean, double std_dev);
+extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed);
 
 //==============================================================================
 //! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
@@ -190,10 +196,12 @@ extern "C" double normal_variate(double mean, double std_dev);
 //! @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 seed A pointer to the pseudorandom seed
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double muir_spectrum(double e0, double m_rat, double kt);
+extern "C" double muir_spectrum(double e0, double m_rat, double kt,
+  uint64_t* seed);
 
 //==============================================================================
 //! Doppler broadens the windowed multipole curvefit.

include/openmc/mgxs.h

@@ -154,8 +154,9 @@ class Mgxs {
     //! @param gin Incoming energy group.
     //! @param dg Sampled delayed group index.
     //! @param gout Sampled outgoing energy group.
+    //! @param seed Pseudorandom seed pointer
     void
-    sample_fission_energy(int gin, int& dg, int& gout);
+    sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed);
 
     //! \brief Samples the outgoing energy and angle from a scatter event.
     //!
@@ -163,8 +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 seed Pseudorandom seed pointer.
     void
-    sample_scatter(int gin, int& gout, double& mu, double& wgt);
+    sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* seed);
 
     //! \brief Calculates cross section quantities needed for tracking.
     //!

include/openmc/particle.h

@@ -13,6 +13,7 @@
 
 #include "openmc/constants.h"
 #include "openmc/position.h"
+#include "openmc/random_lcg.h"
 
 namespace openmc {
 
@@ -217,6 +218,10 @@ public:
   //! create a particle restart HDF5 file
   void write_restart() const;
 
+  //! Gets the pointer to the particle's current PRN seed
+  uint64_t* current_seed() {return seeds_ + stream_;}
+  const uint64_t* current_seed() const {return seeds_ + stream_;}
+
   //==========================================================================
   // Data members
 
@@ -285,6 +290,10 @@ public:
 
   // Track output
   bool write_track_ {false};
+
+  // Current PRNG state
+  uint64_t seeds_[N_STREAMS]; // current seeds
+  int      stream_;           // current RNG stream
 };
 
 } // namespace openmc

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) const;
+    double* mu, int* i_shell, uint64_t* seed) const;
 
-  double rayleigh_scatter(double alpha) const;
+  double rayleigh_scatter(double alpha, uint64_t* seed) const;
 
   void pair_production(double alpha, double* E_electron, double* E_positron,
-    double* mu_electron, double* mu_positron) const;
+    double* mu_electron, double* mu_positron, uint64_t* seed) const;
 
   void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const;
 
@@ -96,14 +96,15 @@ public:
   xt::xtensor<double, 2> dcs_;
 
 private:
-  void compton_doppler(double alpha, double mu, double* E_out, int* i_shell) const;
+  void compton_doppler(double alpha, double mu, double* E_out, int* i_shell,
+                       uint64_t* seed) const;
 };
 
 //==============================================================================
 // Non-member functions
 //==============================================================================
 
-std::pair<double, double> klein_nishina(double alpha);
+std::pair<double, double> klein_nishina(double alpha, uint64_t* seed);
 
 void free_memory_photon();
 

include/openmc/physics.h

@@ -44,7 +44,7 @@ void sample_positron_reaction(Particle* p);
 //!
 //! \param[in] p Particle
 //! \return Index in the data::nuclides vector
-int sample_nuclide(const Particle* p);
+int sample_nuclide(Particle* p);
 
 //! Determine the average total, prompt, and delayed neutrons produced from
 //! fission and creates appropriate bank sites.
@@ -53,13 +53,13 @@ void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 
 int sample_element(Particle* p);
 
-Reaction* sample_fission(int i_nuclide, const Particle* p);
+Reaction* sample_fission(int i_nuclide, Particle* p);
 
-void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_product);
+void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product);
 
 void absorption(Particle* p, int i_nuclide);
 
-void scatter(Particle*, int i_nuclide);
+void scatter(Particle* p, int i_nuclide);
 
 //! Treats the elastic scattering of a neutron with a target.
 void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
@@ -72,15 +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);
+  Direction v_neut, double xs_eff, double kT, uint64_t* 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);
+Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT,
+  uint64_t* seed);
 
-void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site);
+void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in,
+  Particle::Bank* site, uint64_t* seed);
 
 //! handles all reactions with a single secondary neutron (other than fission),
 //! i.e. level scattering, (n,np), (n,na), etc.

include/openmc/plot.h

@@ -14,6 +14,7 @@
 #include "openmc/geometry.h"
 #include "openmc/particle.h"
 #include "openmc/xml_interface.h"
+#include "openmc/random_lcg.h"
 
 namespace openmc {
 
@@ -28,6 +29,9 @@ namespace model {
 extern std::vector<Plot> plots; //!< Plot instance container
 extern std::unordered_map<int, int> plot_map; //!< map of plot ids to index
 
+extern uint64_t plotter_prn_seeds[N_STREAMS]; // Random number seeds used for plotter
+extern int plotter_stream; // Stream index used by the plotter
+
 } // namespace model
 
 //===============================================================================

include/openmc/random_lcg.h

@@ -10,46 +10,66 @@ namespace openmc {
 // Module constants.
 //==============================================================================
 
-extern "C" const int N_STREAMS;
-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 seed Pseudorandom number seed pointer
 //! @return A random number between 0 and 1
 //==============================================================================
 
-extern "C" double prn();
+double prn(uint64_t* 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 seed Pseudorandom number seed
 //! @return A random number between 0 and 1
 //==============================================================================
 
-extern "C" double future_prn(int64_t n);
+double future_prn(int64_t n, uint64_t seed);
 
 //==============================================================================
-//! Set the RNG seed to a unique value 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 id The particle ID
+//! @param offset The offset from the master seed to be used (e.g., for creating
+//! different streams)
+//! @return The initialized seed value
+//==============================================================================
+
+uint64_t init_seed(int64_t id, int offset);
+
+//==============================================================================
+//! 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 seeds Pseudorandom number seed array
 //! @param id The particle ID
 //==============================================================================
 
-extern "C" void set_particle_seed(int64_t id);
+void init_particle_seeds(int64_t id, uint64_t* seeds);
 
 //==============================================================================
-//! Advance the random number seed 'n' times from the current seed.
+//! 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 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);
+void advance_prn_seed(int64_t n, uint64_t* seed);
 
 //==============================================================================
 //! Advance a random number seed 'n' times.
@@ -63,18 +83,6 @@ extern "C" void advance_prn_seed(int64_t n);
 
 uint64_t future_seed(uint64_t n, uint64_t seed);
 
-//==============================================================================
-//! Switch the RNG to a different stream of random numbers.
-//!
-//! If random numbers are needed in routines not used directly for tracking
-//! (e.g. physics), this allows the numbers to be generated without affecting
-//! reproducibility of the physics.
-//! @param n The RNG stream to switch to. Use the constants such as
-//! `STREAM_TRACKING` and `STREAM_TALLIES` for this argument.
-//==============================================================================
-
-extern "C" void prn_set_stream(int n);
-
 //==============================================================================
 //                               API FUNCTIONS
 //==============================================================================

include/openmc/reaction_product.h

@@ -42,7 +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
-  void sample(double E_in, double& E_out, double& mu) const;
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
   Particle::Type particle_; //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode

include/openmc/scattdata.h

@@ -65,8 +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 seed Pseudorandom number seed pointer
     virtual void
-    sample(int gin, int& gout, double& mu, double& wgt) = 0;
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed) = 0;
 
     //! \brief Initializes the ScattData object from a given scatter and
     //!   multiplicity matrix.
@@ -109,8 +110,9 @@ class ScattData {
     //! @param gin Incoming energy group.
     //! @param gout Sampled outgoing energy group.
     //! @param i_gout Sampled outgoing energy group index.
+    //! @param seed Pseudorandom number seed pointer
     void
-    sample_energy(int gin, int& gout, int& i_gout);
+    sample_energy(int gin, int& gout, int& i_gout, uint64_t* seed);
 
     //! \brief Provides a cross section value given certain parameters
     //!
@@ -161,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);
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed);
 
     size_t
     get_order() {return dist[0][0].size() - 1;};
@@ -197,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);
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed);
 
     size_t
     get_order() {return dist[0][0].size();};
@@ -238,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);
+    sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed);
 
     size_t
     get_order() {return dist[0][0].size();};

include/openmc/secondary_correlated.h

@@ -38,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 
   // energy property
   std::vector<double>& energy() { return energy_; }

include/openmc/secondary_kalbach.h

@@ -29,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 private:
   //! Outgoing energy/angle at a single incoming energy
   struct KMTable {

include/openmc/secondary_nbody.h

@@ -24,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 private:
   int n_bodies_; //!< Number of particles distributed
   double mass_ratio_; //!< Total mass of particles [neutron mass]

include/openmc/secondary_thermal.h

@@ -31,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
 };
@@ -51,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 private:
   double debye_waller_;
 };
@@ -72,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* 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
@@ -94,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 private:
   const std::vector<double>& energy_; //!< Incident energies
   xt::xtensor<double, 2> energy_out_; //!< Outgoing energies for each incident energy
@@ -117,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom number seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 private:
   //! Secondary energy/angle distribution
   struct DistEnergySab {

include/openmc/secondary_uncorrelated.h

@@ -29,7 +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
-  void sample(double E_in, double& E_out, double& mu) const override;
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample(double E_in, double& E_out, double& mu,
+    uint64_t* seed) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }

include/openmc/source.h

@@ -38,8 +38,9 @@ public:
   explicit SourceDistribution(pugi::xml_node node);
 
   //! Sample from the external source distribution
+  //! \param[inout] seed Pseudorandom seed pointer
   //! \return Sampled site
-  Particle::Bank sample() const;
+  Particle::Bank sample(uint64_t* seed) const;
 
   // Properties
   double strength() const { return strength_; }
@@ -60,8 +61,9 @@ extern "C" void initialize_source();
 
 //! Sample a site from all external source distributions in proportion to their
 //! source strength
+//! \param[inout] seed Pseudorandom seed pointer
 //! \return Sampled source site
-Particle::Bank sample_external_source();
+Particle::Bank sample_external_source(uint64_t* seed);
 
 //! Fill source bank at end of generation for fixed source simulations
 void fill_source_bank_fixedsource();

include/openmc/surface.h

@@ -116,7 +116,8 @@ public:
   //! \return Outgoing direction of the ray
   virtual Direction reflect(Position r, Direction u) const;
 
-  virtual Direction diffuse_reflect(Position r, Direction u) const;
+  virtual Direction diffuse_reflect(Position r, Direction u,
+    uint64_t* seed) const;
 
   //! Evaluate the equation describing the surface.
   //!

include/openmc/thermal.h

@@ -64,8 +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] seed Pseudorandom seed pointer
   void sample(const NuclideMicroXS& micro_xs, double E_in,
-              double* E_out, double* mu);
+              double* E_out, double* mu, uint64_t* seed);
 private:
   struct Reaction {
     // Default constructor
@@ -100,8 +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] seed Pseudorandom seed pointer
   void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic,
-                    double* inelastic) const;
+                    double* inelastic, uint64_t* seed) const;
 
   //! Determine whether table applies to a particular nuclide
   //!