diff --git a/CMakeLists.txt b/CMakeLists.txt
index 87b8789d101..d3119fb8758 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -372,6 +372,7 @@ list(APPEND libopenmc_SOURCES
   src/particle.cpp
   src/particle_data.cpp
   src/particle_restart.cpp
+  src/particle_type.cpp
   src/photon.cpp
   src/physics.cpp
   src/physics_common.cpp
diff --git a/docs/source/io_formats/collision_track.rst b/docs/source/io_formats/collision_track.rst
index a3645975444..8e1e00ffb8e 100644
--- a/docs/source/io_formats/collision_track.rst
+++ b/docs/source/io_formats/collision_track.rst
@@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested
 (``collision_track.N.h5``) or when the run is performed in parallel. The file
 contains the information needed to reconstruct each recorded collision.
 
-The current revision of the collision track file format is 1.0.
+The current revision of the collision track file format is 1.1.
 
 **/**
 
@@ -37,9 +37,9 @@ The current revision of the collision track file format is 1.0.
               - ``material_id`` (*int*) -- ID of the material containing the collision site.
               - ``universe_id`` (*int*) -- ID of the universe containing the collision site.
               - ``n_collision`` (*int*) -- Collision counter for the particle history.
-              - ``particle`` (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, 3=positron).
-              - ``parent_id`` (*int64*) -- Unique ID of the parent particle.
-              - ``progeny_id`` (*int64*) -- Progeny ID of the particle.
+              - ``particle`` (*int32_t*) -- Particle type (PDG number).
+              - ``parent_id`` (*int64_t*) -- Unique ID of the parent particle.
+              - ``progeny_id`` (*int64_t*) -- Progeny ID of the particle.
 
 In an MPI run, OpenMC writes the combined dataset by gathering collision-track
 entries from all ranks before flushing them to disk, so the final file appears
diff --git a/docs/source/io_formats/particle_restart.rst b/docs/source/io_formats/particle_restart.rst
index 2734f0470ef..e0fe76b9d81 100644
--- a/docs/source/io_formats/particle_restart.rst
+++ b/docs/source/io_formats/particle_restart.rst
@@ -4,7 +4,7 @@
 Particle Restart File Format
 ============================
 
-The current version of the particle restart file format is 2.0.
+The current version of the particle restart file format is 2.1.
 
 **/**
 
@@ -26,8 +26,7 @@ The current version of the particle restart file format is 2.0.
            - **run_mode** (*char[]*) -- Run mode used, either 'fixed source',
              'eigenvalue', or 'particle restart'.
            - **id** (*int8_t*) -- Unique identifier of the particle.
-           - **type** (*int*) -- Particle type (0=neutron, 1=photon, 2=electron,
-             3=positron)
+           - **type** (*int32_t*) -- Particle type (PDG number)
            - **weight** (*double*) -- Weight of the particle.
            - **energy** (*double*) -- Energy of the particle in eV for
              continuous-energy mode, or the energy group of the particle for
diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst
index 8fc26ae070b..ed7c6273aaa 100644
--- a/docs/source/io_formats/settings.rst
+++ b/docs/source/io_formats/settings.rst
@@ -721,7 +721,10 @@ attributes/sub-elements:
     is present.
 
   :particle:
-    The source particle type, either ``neutron`` or ``photon``.
+    The source particle type, specified as a PDG number or a string alias (e.g.,
+    ``neutron``/``n``, ``photon``/``gamma``, ``electron``, ``positron``,
+    ``proton``/``p``, ``deuteron``/``d``, ``triton``/``t``, ``alpha``, or GNDS
+    nuclide names like ``Fe57``).
 
     *Default*: neutron
 
@@ -1537,7 +1540,8 @@ sub-elements/attributes:
     *Default*: None
 
   :particle_type:
-    The particle that the weight windows will apply to (e.g., 'neutron')
+    The particle that the weight windows will apply to, specified as a PDG
+    code or string (e.g., ``neutron``).
 
     *Default*: 'neutron'
 
@@ -1597,7 +1601,8 @@ mesh-based weight windows.
     *Default*: None
 
   :particle_type:
-    The particle that the weight windows will apply to (e.g., 'neutron')
+    The particle that the weight windows will apply to, specified as a PDG
+    code or string (e.g., ``neutron``).
 
     *Default*: neutron
 
diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst
index 4ce2229ed8f..2e81e0a3f2f 100644
--- a/docs/source/io_formats/source.rst
+++ b/docs/source/io_formats/source.rst
@@ -15,6 +15,8 @@ following the same format.
 **/**
 
 :Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
+             - **version** (*int[2]*) -- Major and minor version of the source
+               file format.
 
 :Datasets:
 
@@ -22,5 +24,5 @@ following the same format.
              particle. The compound type has fields ``r``, ``u``, ``E``,
              ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``,
              which represent the position, direction, energy, time, weight,
-             delayed group, surface ID, and particle type (0=neutron, 1=photon,
-             2=electron, 3=positron), respectively.
+             delayed group, surface ID, and particle type (PDG number),
+             respectively.
diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst
index 2309643dc89..7d7765b849d 100644
--- a/docs/source/io_formats/statepoint.rst
+++ b/docs/source/io_formats/statepoint.rst
@@ -4,7 +4,7 @@
 State Point File Format
 =======================
 
-The current version of the statepoint file format is 18.1.
+The current version of the statepoint file format is 18.2.
 
 **/**
 
@@ -56,8 +56,8 @@ The current version of the statepoint file format is 18.1.
              ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and
              ``particle``, which represent the position, direction, energy,
              time, weight, delayed group, surface ID, and particle type
-             (0=neutron, 1=photon, 2=electron, 3=positron), respectively. Only
-             present when `run_mode` is 'eigenvalue'.
+             (PDG number), respectively. Only present when `run_mode` is
+             'eigenvalue'.
 
 **/tallies/**
 
diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst
index 9f29a949acb..0ba0e061fd5 100644
--- a/docs/source/io_formats/tallies.rst
+++ b/docs/source/io_formats/tallies.rst
@@ -318,8 +318,8 @@ should be set to:
   they use ``energy`` and ``y``.
 
 :particle:
-  A list of integers indicating the type of particles to tally ('neutron' = 1,
-  'photon' = 2, 'electron' = 3, 'positron' = 4).
+  A list of particle identifiers to tally, specified as strings (e.g.,
+  ``neutron``, ``photon``, ``He4``) or as integer PDG numbers.
 
 ------------------
 ``<mesh>`` Element
diff --git a/docs/source/io_formats/track.rst b/docs/source/io_formats/track.rst
index a97d75e58e6..9c05ac4c370 100644
--- a/docs/source/io_formats/track.rst
+++ b/docs/source/io_formats/track.rst
@@ -4,7 +4,7 @@
 Track File Format
 =================
 
-The current revision of the particle track file format is 3.0.
+The current revision of the particle track file format is 3.1.
 
 **/**
 
@@ -32,6 +32,5 @@ The current revision of the particle track file format is 3.0.
                             the array for each primary/secondary particle. The
                             last offset should match the total size of the
                             array.
-                          - **particles** (*int[]*) -- Particle type for each
-                            primary/secondary particle (0=neutron, 1=photon,
-                            2=electron, 3=positron).
+                          - **particles** (*int32_t[]*) -- Particle type for
+                            each primary/secondary particle (PDG number).
diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst
index b9c693ba1cd..5a04fedd70a 100644
--- a/docs/source/usersguide/settings.rst
+++ b/docs/source/usersguide/settings.rst
@@ -400,7 +400,7 @@ below.
     {
       openmc::SourceSite particle;
       // weight
-      particle.particle = openmc::ParticleType::neutron;
+      particle.particle = openmc::ParticleType::neutron();
       particle.wgt = 1.0;
       // position
       double angle = 2.0 * M_PI * openmc::prn(seed);
@@ -477,7 +477,7 @@ parameters to the source class when it is created:
     {
       openmc::SourceSite particle;
       // weight
-      particle.particle = openmc::ParticleType::neutron;
+      particle.particle = openmc::ParticleType::neutron();
       particle.wgt = 1.0;
       // position
       particle.r.x = 0.0;
diff --git a/examples/custom_source/source_ring.cpp b/examples/custom_source/source_ring.cpp
index b87afbd1760..5ab531a4a7d 100644
--- a/examples/custom_source/source_ring.cpp
+++ b/examples/custom_source/source_ring.cpp
@@ -1,17 +1,16 @@
-#include <cmath> // for M_PI
+#include <cmath>  // for M_PI
 #include <memory> // for unique_ptr
 
+#include "openmc/particle.h"
 #include "openmc/random_lcg.h"
 #include "openmc/source.h"
-#include "openmc/particle.h"
 
-class RingSource : public openmc::Source
-{
+class RingSource : public openmc::Source {
   openmc::SourceSite sample(uint64_t* seed) const
   {
     openmc::SourceSite particle;
     // particle type
-    particle.particle = openmc::ParticleType::neutron;
+    particle.particle = openmc::ParticleType::neutron();
     // position
     double angle = 2.0 * M_PI * openmc::prn(seed);
     double radius = 3.0;
@@ -25,10 +24,11 @@ class RingSource : public openmc::Source
   }
 };
 
-// A function to create a unique pointer to an instance of this class when generated
-// via a plugin call using dlopen/dlsym.
-// You must have external C linkage here otherwise dlopen will not find the file
-extern "C" std::unique_ptr<RingSource> openmc_create_source(std::string parameters)
+// A function to create a unique pointer to an instance of this class when
+// generated via a plugin call using dlopen/dlsym. You must have external C
+// linkage here otherwise dlopen will not find the file
+extern "C" std::unique_ptr<RingSource> openmc_create_source(
+  std::string parameters)
 {
   return std::make_unique<RingSource>();
 }
diff --git a/examples/parameterized_custom_source/parameterized_source_ring.cpp b/examples/parameterized_custom_source/parameterized_source_ring.cpp
index e70dc8bb190..a983414f1c5 100644
--- a/examples/parameterized_custom_source/parameterized_source_ring.cpp
+++ b/examples/parameterized_custom_source/parameterized_source_ring.cpp
@@ -1,63 +1,65 @@
-#include <cmath> // for M_PI
+#include <cmath>  // for M_PI
 #include <memory> // for unique_ptr
 #include <unordered_map>
 
+#include "openmc/particle.h"
 #include "openmc/random_lcg.h"
 #include "openmc/source.h"
-#include "openmc/particle.h"
 
 class RingSource : public openmc::Source {
-  public:
-    RingSource(double radius, double energy) : radius_(radius), energy_(energy) { }
-
-    // Defines a function that can create a unique pointer to a new instance of this class
-    // by extracting the parameters from the provided string.
-    static std::unique_ptr<RingSource> from_string(std::string parameters)
-    {
-      std::unordered_map<std::string, std::string> parameter_mapping;
-
-      std::stringstream ss(parameters);
-      std::string parameter;
-      while (std::getline(ss, parameter, ',')) {
-        parameter.erase(0, parameter.find_first_not_of(' '));
-        std::string key = parameter.substr(0, parameter.find_first_of('='));
-        std::string value = parameter.substr(parameter.find_first_of('=') + 1, parameter.length());
-        parameter_mapping[key] = value;
-      }
-
-      double radius = std::stod(parameter_mapping["radius"]);
-      double energy = std::stod(parameter_mapping["energy"]);
-      return std::make_unique<RingSource>(radius, energy);
-    }
+public:
+  RingSource(double radius, double energy) : radius_(radius), energy_(energy) {}
+
+  // Defines a function that can create a unique pointer to a new instance of
+  // this class by extracting the parameters from the provided string.
+  static std::unique_ptr<RingSource> from_string(std::string parameters)
+  {
+    std::unordered_map<std::string, std::string> parameter_mapping;
 
-    // Samples from an instance of this class.
-    openmc::SourceSite sample(uint64_t* seed) const
-    {
-      openmc::SourceSite particle;
-      // particle type
-      particle.particle = openmc::ParticleType::neutron;
-      // position
-      double angle = 2.0 * M_PI * openmc::prn(seed);
-      double radius = this->radius_;
-      particle.r.x = radius * std::cos(angle);
-      particle.r.y = radius * std::sin(angle);
-      particle.r.z = 0.0;
-      // angle
-      particle.u = {1.0, 0.0, 0.0};
-      particle.E = this->energy_;
-
-      return particle;
+    std::stringstream ss(parameters);
+    std::string parameter;
+    while (std::getline(ss, parameter, ',')) {
+      parameter.erase(0, parameter.find_first_not_of(' '));
+      std::string key = parameter.substr(0, parameter.find_first_of('='));
+      std::string value =
+        parameter.substr(parameter.find_first_of('=') + 1, parameter.length());
+      parameter_mapping[key] = value;
     }
 
-  private:
-    double radius_;
-    double energy_;
+    double radius = std::stod(parameter_mapping["radius"]);
+    double energy = std::stod(parameter_mapping["energy"]);
+    return std::make_unique<RingSource>(radius, energy);
+  }
+
+  // Samples from an instance of this class.
+  openmc::SourceSite sample(uint64_t* seed) const
+  {
+    openmc::SourceSite particle;
+    // particle type
+    particle.particle = openmc::ParticleType::neutron();
+    // position
+    double angle = 2.0 * M_PI * openmc::prn(seed);
+    double radius = this->radius_;
+    particle.r.x = radius * std::cos(angle);
+    particle.r.y = radius * std::sin(angle);
+    particle.r.z = 0.0;
+    // angle
+    particle.u = {1.0, 0.0, 0.0};
+    particle.E = this->energy_;
+
+    return particle;
+  }
+
+private:
+  double radius_;
+  double energy_;
 };
 
-// A function to create a unique pointer to an instance of this class when generated
-// via a plugin call using dlopen/dlsym.
-// You must have external C linkage here otherwise dlopen will not find the file
-extern "C" std::unique_ptr<RingSource> openmc_create_source(std::string parameters)
+// A function to create a unique pointer to an instance of this class when
+// generated via a plugin call using dlopen/dlsym. You must have external C
+// linkage here otherwise dlopen will not find the file
+extern "C" std::unique_ptr<RingSource> openmc_create_source(
+  std::string parameters)
 {
   return RingSource::from_string(parameters);
 }
diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index d8041ef4149..3ce7407c1a7 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -201,8 +201,8 @@ int openmc_weight_windows_set_energy_bounds(
   int32_t index, double* e_bounds, size_t e_bounds_size);
 int openmc_weight_windows_get_energy_bounds(
   int32_t index, const double** e_bounds, size_t* e_bounds_size);
-int openmc_weight_windows_set_particle(int32_t index, int particle);
-int openmc_weight_windows_get_particle(int32_t index, int* particle);
+int openmc_weight_windows_set_particle(int32_t index, int32_t particle);
+int openmc_weight_windows_get_particle(int32_t index, int32_t* particle);
 int openmc_weight_windows_get_bounds(int32_t index, const double** lower_bounds,
   const double** upper_bounds, size_t* size);
 int openmc_weight_windows_set_bounds(int32_t index, const double* lower_bounds,
@@ -227,7 +227,7 @@ int openmc_zernike_filter_set_order(int32_t index, int order);
 int openmc_zernike_filter_set_params(
   int32_t index, const double* x, const double* y, const double* r);
 
-int openmc_particle_filter_get_bins(int32_t idx, int bins[]);
+int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]);
 
 //! Sets the mesh and energy grid for CMFD reweight
 //! \param[in] meshtyally_id id of CMFD Mesh Tally
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index bba260c3a4b..980cff447ef 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -25,16 +25,16 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
 constexpr int HDF5_VERSION[] {3, 0};
 
 // Version numbers for binary files
-constexpr array<int, 2> VERSION_STATEPOINT {18, 1};
-constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
-constexpr array<int, 2> VERSION_TRACK {3, 0};
+constexpr array<int, 2> VERSION_STATEPOINT {18, 2};
+constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 1};
+constexpr array<int, 2> VERSION_TRACK {3, 1};
 constexpr array<int, 2> VERSION_SUMMARY {6, 1};
 constexpr array<int, 2> VERSION_VOLUME {1, 0};
 constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
 constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
 constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
-constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 0};
+constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 1};
 
 // ============================================================================
 // ADJUSTABLE PARAMETERS
diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h
index 60b88a153bc..58d83393963 100644
--- a/include/openmc/nuclide.h
+++ b/include/openmc/nuclide.h
@@ -163,7 +163,7 @@ bool multipole_in_range(const Nuclide& nuc, double E);
 namespace data {
 
 // Minimum/maximum transport energy for each particle type. Order corresponds to
-// that of the ParticleType enum
+// transport_index() for supported transport particles.
 extern array<double, 4> energy_min;
 extern array<double, 4> energy_max;
 
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 0f37719b94f..2f6e6196bf1 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -126,10 +126,6 @@ class Particle : public ParticleData {
 //! Functions
 //============================================================================
 
-std::string particle_type_to_str(ParticleType type);
-
-ParticleType str_to_particle_type(std::string str);
-
 void add_surf_source_to_bank(Particle& p, const Surface& surf);
 
 } // namespace openmc
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index fdacfa765b3..75166b36af5 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -3,6 +3,7 @@
 
 #include "openmc/array.h"
 #include "openmc/constants.h"
+#include "openmc/particle_type.h"
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
 #include "openmc/tallies/filter_match.h"
@@ -30,9 +31,6 @@ constexpr double CACHE_INVALID {-1.0};
 //==========================================================================
 // Aliases and type definitions
 
-//! Particle types
-enum class ParticleType { neutron, photon, electron, positron };
-
 //! Saved ("banked") state of a particle
 //! NOTE: This structure's MPI type is built in initialize_mpi() of
 //! initialize.cpp. Any changes made to the struct here must also be
@@ -496,7 +494,7 @@ class ParticleData : public GeometryState {
   MacroXS macro_xs_;
   CacheDataMG mg_xs_cache_;
 
-  ParticleType type_ {ParticleType::neutron};
+  ParticleType type_;
 
   double E_;
   double E_last_;
diff --git a/include/openmc/particle_type.h b/include/openmc/particle_type.h
new file mode 100644
index 00000000000..4cd48d047bd
--- /dev/null
+++ b/include/openmc/particle_type.h
@@ -0,0 +1,172 @@
+//==============================================================================
+// ParticleType class definition
+//==============================================================================
+
+#ifndef OPENMC_PARTICLE_TYPE_H
+#define OPENMC_PARTICLE_TYPE_H
+
+#include <cstddef>
+#include <cstdint>
+#include <string>
+#include <string_view>
+#include <type_traits>
+
+#include "openmc/constants.h"
+
+namespace openmc {
+
+//------------------------------------------------------------------------------
+// PDG constants (canonical particle identity as simple integers)
+//------------------------------------------------------------------------------
+
+inline constexpr int32_t PDG_NEUTRON = 2112;
+inline constexpr int32_t PDG_PHOTON = 22;
+inline constexpr int32_t PDG_ELECTRON = 11;
+inline constexpr int32_t PDG_POSITRON = -11;
+inline constexpr int32_t PDG_PROTON = 2212;
+inline constexpr int32_t PDG_DEUTERON = 1000010020;
+inline constexpr int32_t PDG_TRITON = 1000010030;
+inline constexpr int32_t PDG_ALPHA = 1000020040;
+
+//------------------------------------------------------------------------------
+// ParticleType class (standard-layout, trivially copyable)
+//------------------------------------------------------------------------------
+
+class ParticleType {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors
+
+  // Default constructor: defaults to neutron
+  constexpr ParticleType() : pdg_number_(PDG_NEUTRON) {}
+
+  // Constructor from PDG number
+  constexpr explicit ParticleType(int32_t pdg_number) : pdg_number_(pdg_number)
+  {}
+
+  // Constructor from particle name string (e.g., "neutron", "photon", "Fe56")
+  explicit ParticleType(std::string_view str);
+
+  // Constructor from Z, A, and metastable state for nuclear particles
+  constexpr ParticleType(int Z, int A, int m = 0)
+    : pdg_number_(1000000000 + Z * 10000 + A * 10 + m)
+  {}
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  // Accessor for the underlying PDG number
+  constexpr int32_t pdg_number() const { return pdg_number_; }
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  // Convert to string representation
+  std::string str() const;
+
+  // Check if this represents a nucleus (vs elementary particle)
+  constexpr bool is_nucleus() const
+  {
+    // PDG nuclear codes are >= 1000000000 (100ZZZAAAI format)
+    return pdg_number_ >= 1000000000;
+  }
+
+  // Get transport index (0-3 for transportable particles, C_NONE otherwise)
+  constexpr int transport_index() const;
+
+  // Check if this is a neutron
+  constexpr bool is_neutron() const { return pdg_number_ == PDG_NEUTRON; }
+
+  // Check if this is a photon
+  constexpr bool is_photon() const { return pdg_number_ == PDG_PHOTON; }
+
+  constexpr bool is_transportable()
+  {
+    return this->transport_index() != C_NONE;
+  }
+
+  //----------------------------------------------------------------------------
+  // Static factory methods
+
+  static constexpr ParticleType neutron() { return ParticleType {PDG_NEUTRON}; }
+  static constexpr ParticleType photon() { return ParticleType {PDG_PHOTON}; }
+  static constexpr ParticleType electron()
+  {
+    return ParticleType {PDG_ELECTRON};
+  }
+  static constexpr ParticleType positron()
+  {
+    return ParticleType {PDG_POSITRON};
+  }
+  static constexpr ParticleType proton() { return ParticleType {PDG_PROTON}; }
+  static constexpr ParticleType deuteron()
+  {
+    return ParticleType {PDG_DEUTERON};
+  }
+  static constexpr ParticleType triton() { return ParticleType {PDG_TRITON}; }
+  static constexpr ParticleType alpha() { return ParticleType {PDG_ALPHA}; }
+
+private:
+  int32_t pdg_number_;
+};
+
+//------------------------------------------------------------------------------
+// Static assertions to ensure standard-layout and trivially copyable
+//------------------------------------------------------------------------------
+
+static_assert(std::is_standard_layout_v<ParticleType>,
+  "ParticleType must be standard-layout");
+static_assert(std::is_trivially_copyable_v<ParticleType>,
+  "ParticleType must be trivially copyable");
+static_assert(sizeof(ParticleType) == sizeof(int32_t),
+  "ParticleType must be same size as int32_t");
+
+//------------------------------------------------------------------------------
+// Comparison operators (free functions for symmetry)
+//------------------------------------------------------------------------------
+
+constexpr bool operator==(ParticleType lhs, ParticleType rhs)
+{
+  return lhs.pdg_number() == rhs.pdg_number();
+}
+
+constexpr bool operator!=(ParticleType lhs, ParticleType rhs)
+{
+  return lhs.pdg_number() != rhs.pdg_number();
+}
+
+constexpr bool operator<(ParticleType lhs, ParticleType rhs)
+{
+  return lhs.pdg_number() < rhs.pdg_number();
+}
+
+//------------------------------------------------------------------------------
+// ParticleType member function implementations (inline)
+//------------------------------------------------------------------------------
+
+constexpr int ParticleType::transport_index() const
+{
+  switch (pdg_number_) {
+  case PDG_NEUTRON:
+    return 0;
+  case PDG_PHOTON:
+    return 1;
+  case PDG_ELECTRON:
+    return 2;
+  case PDG_POSITRON:
+    return 3;
+  default:
+    return C_NONE;
+  }
+}
+
+//------------------------------------------------------------------------------
+// Legacy conversion helpers
+//------------------------------------------------------------------------------
+
+// Legacy enum code (0..3) to ParticleType conversion
+ParticleType legacy_particle_index_to_type(int code);
+
+} // namespace openmc
+
+#endif // OPENMC_PARTICLE_TYPE_H
diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h
index 4fbbc1b626a..9a8eab7d913 100644
--- a/include/openmc/reaction_product.h
+++ b/include/openmc/reaction_product.h
@@ -10,7 +10,7 @@
 #include "openmc/chain.h"
 #include "openmc/endf.h"
 #include "openmc/memory.h" // for unique_ptr
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/vector.h" // for vector
 
 namespace openmc {
diff --git a/include/openmc/source.h b/include/openmc/source.h
index 36c821fc07c..2f32aa2a058 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -12,7 +12,7 @@
 #include "openmc/distribution_multi.h"
 #include "openmc/distribution_spatial.h"
 #include "openmc/memory.h"
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -148,11 +148,11 @@ class IndependentSource : public Source {
 
 private:
   // Data members
-  ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted
-  UPtrSpace space_;                               //!< Spatial distribution
-  UPtrAngle angle_;                               //!< Angular distribution
-  UPtrDist energy_;                               //!< Energy distribution
-  UPtrDist time_;                                 //!< Time distribution
+  ParticleType particle_; //!< Type of particle emitted
+  UPtrSpace space_;       //!< Spatial distribution
+  UPtrAngle angle_;       //!< Angular distribution
+  UPtrDist energy_;       //!< Energy distribution
+  UPtrDist time_;         //!< Time distribution
 };
 
 //==============================================================================
diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h
index 863a6d282fe..9932a6037e7 100644
--- a/include/openmc/tallies/filter_particle.h
+++ b/include/openmc/tallies/filter_particle.h
@@ -1,7 +1,7 @@
 #ifndef OPENMC_TALLIES_FILTER_PARTICLE_H
 #define OPENMC_TALLIES_FILTER_PARTICLE_H
 
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/vector.h"
diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h
index 7638155228e..0d7435ca2a5 100644
--- a/include/openmc/weight_windows.h
+++ b/include/openmc/weight_windows.h
@@ -10,7 +10,7 @@
 #include "openmc/constants.h"
 #include "openmc/memory.h"
 #include "openmc/mesh.h"
-#include "openmc/particle.h"
+#include "openmc/particle_type.h"
 #include "openmc/span.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/vector.h"
@@ -193,10 +193,9 @@ class WeightWindows {
 private:
   //----------------------------------------------------------------------------
   // Data members
-  int32_t id_;    //!< Unique ID
-  int64_t index_; //!< Index into weight windows vector
-  ParticleType particle_type_ {
-    ParticleType::neutron};      //!< Particle type to apply weight windows to
+  int32_t id_;                   //!< Unique ID
+  int64_t index_;                //!< Index into weight windows vector
+  ParticleType particle_type_;   //!< Particle type to apply weight windows to
   vector<double> energy_bounds_; //!< Energy boundaries [eV]
   xt::xtensor<double, 2> lower_ww_; //!< Lower weight window bounds (shape:
                                     //!< energy_bins, mesh_bins (k, j, i))
diff --git a/openmc/filter.py b/openmc/filter.py
index e7c54f2f0ee..53ec93e21d1 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -35,7 +35,6 @@
     'z-min out', 'z-min in', 'z-max out', 'z-max in'
 )
 
-_PARTICLES = {'neutron', 'photon', 'electron', 'positron'}
 
 
 class FilterMeta(ABCMeta):
@@ -735,9 +734,8 @@ class ParticleFilter(Filter):
 
     Parameters
     ----------
-    bins : str, or sequence of str
-        The particles to tally represented as strings ('neutron', 'photon',
-        'electron', 'positron').
+    bins : str, int, openmc.ParticleType, or sequence
+        The particle types to tally represented as names, PDG numbers, or types.
     filter_id : int
         Unique identifier for the filter
 
@@ -763,11 +761,16 @@ def __eq__(self, other):
 
     @Filter.bins.setter
     def bins(self, bins):
-        cv.check_type('bins', bins, Sequence, str)
+        if isinstance(bins, (str, Integral, openmc.ParticleType)):
+            bins = [bins]
+        else:
+            cv.check_type('bins', bins, Sequence,
+                          (str, Integral, openmc.ParticleType))
         bins = np.atleast_1d(bins)
-        for edge in bins:
-            cv.check_value('filter bin', edge, _PARTICLES)
-        self._bins = bins
+        normalized = []
+        for entry in bins:
+            normalized.append(str(openmc.ParticleType(entry)))
+        self._bins = np.array(normalized, dtype=str)
 
     @classmethod
     def from_hdf5(cls, group, **kwargs):
diff --git a/openmc/lib/core.py b/openmc/lib/core.py
index 02c7784d1cc..79310bb37f9 100644
--- a/openmc/lib/core.py
+++ b/openmc/lib/core.py
@@ -28,7 +28,7 @@ class _SourceSite(Structure):
                 ('wgt', c_double),
                 ('delayed_group', c_int),
                 ('surf_id', c_int),
-                ('particle', c_int),
+                ('particle', c_int32),
                 ('parent_nuclide', c_int),
                 ('parent_id', c_int64),
                 ('progeny_id', c_int64)]
diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py
index 5cf496107b3..dc81c3487b9 100644
--- a/openmc/lib/filter.py
+++ b/openmc/lib/filter.py
@@ -132,6 +132,9 @@
 _dll.openmc_new_filter.argtypes = [c_char_p, POINTER(c_int32)]
 _dll.openmc_new_filter.restype = c_int
 _dll.openmc_new_filter.errcheck = _error_handler
+_dll.openmc_particle_filter_get_bins.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_particle_filter_get_bins.restype = c_int
+_dll.openmc_particle_filter_get_bins.errcheck = _error_handler
 _dll.openmc_spatial_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)]
 _dll.openmc_spatial_legendre_filter_get_order.restype = c_int
 _dll.openmc_spatial_legendre_filter_get_order.errcheck = _error_handler
@@ -402,8 +405,8 @@ class MeshFilter(Filter):
     translation : Iterable of float
         3-D coordinates of the translation vector
     rotation : Iterable of float
-        The rotation matrix or angles of the filter mesh. This can either be 
-        a fully specified 3 x 3 rotation matrix or an Iterable of length 3 
+        The rotation matrix or angles of the filter mesh. This can either be
+        a fully specified 3 x 3 rotation matrix or an Iterable of length 3
         with the angles in degrees about the x, y, and z axes, respectively.
 
     """
@@ -454,7 +457,7 @@ def rotation(self):
         else:
             raise ValueError(
                 f'Invalid size of rotation matrix: {rot_size}')
-        
+
     @rotation.setter
     def rotation(self, rotation_data):
         flat_rotation = np.asarray(rotation_data, dtype=float).flatten()
@@ -598,9 +601,9 @@ class ParticleFilter(Filter):
 
     @property
     def bins(self):
-        particle_i = np.zeros((self.n_bins,), dtype=c_int)
+        particle_i = np.zeros((self.n_bins,), dtype=c_int32)
         _dll.openmc_particle_filter_get_bins(
-            self._index, particle_i.ctypes.data_as(POINTER(c_int)))
+            self._index, particle_i.ctypes.data_as(POINTER(c_int32)))
         return [ParticleType(i) for i in particle_i]
 
 
diff --git a/openmc/lib/weight_windows.py b/openmc/lib/weight_windows.py
index ed442d33ff4..2b26d3b55f5 100644
--- a/openmc/lib/weight_windows.py
+++ b/openmc/lib/weight_windows.py
@@ -53,11 +53,11 @@
 _dll.openmc_weight_windows_get_energy_bounds.restype = c_int
 _dll.openmc_weight_windows_get_energy_bounds.errcheck = _error_handler
 
-_dll.openmc_weight_windows_set_particle.argtypes = [c_int32, c_int]
+_dll.openmc_weight_windows_set_particle.argtypes = [c_int32, c_int32]
 _dll.openmc_weight_windows_set_particle.restype = c_int
 _dll.openmc_weight_windows_set_particle.errcheck = _error_handler
 
-_dll.openmc_weight_windows_get_particle.argtypes = [c_int32, POINTER(c_int)]
+_dll.openmc_weight_windows_get_particle.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_weight_windows_get_particle.restype = c_int
 _dll.openmc_weight_windows_get_particle.errcheck = _error_handler
 
@@ -201,16 +201,13 @@ def energy_bounds(self, e_bounds):
 
     @property
     def particle(self):
-        val = c_int()
+        val = c_int32()
         _dll.openmc_weight_windows_get_particle(self._index, val)
         return ParticleType(val.value)
 
     @particle.setter
     def particle(self, p):
-        if isinstance(p, str):
-            p = ParticleType.from_string(p)
-        else:
-            p = ParticleType(p)
+        p = ParticleType(p)
         _dll.openmc_weight_windows_set_particle(self._index, int(p))
 
     @property
@@ -304,10 +301,10 @@ def from_tally(cls, tally, particle=ParticleType.NEUTRON):
         ----------
         tally : openmc.lib.Tally
             The tally used to create the WeightWindows instance.
-        particle : openmc.ParticleType or str, optional
+        particle : openmc.ParticleType or str or int, optional
             The particle type to use for the WeightWindows instance. Should be
-            specified as an instance of ParticleType or as a string with a value of
-            'neutron' or 'photon'.
+            specified as an instance of ParticleType, a PDG number, or as a
+            name.
 
         Returns
         -------
@@ -317,7 +314,8 @@ def from_tally(cls, tally, particle=ParticleType.NEUTRON):
         Raises
         ------
         ValueError
-            If the particle parameter is not an instance of ParticleType or a string.
+            If the particle parameter is not an instance of ParticleType, a string,
+            or an integer PDG number.
         ValueError
             If the particle parameter is not a valid particle type (i.e., not 'neutron'
             or 'photon').
@@ -328,12 +326,13 @@ def from_tally(cls, tally, particle=ParticleType.NEUTRON):
             If the tally does not have a MeshFilter.
         """
         # do some checks on particle value
-        if not isinstance(particle, (ParticleType, str)):
-            raise ValueError(f"Parameter 'particle' must be {ParticleType} or one of ('neutron', 'photon').")
+        if not isinstance(particle, (ParticleType, str, int)):
+            raise ValueError(
+                f"Parameter 'particle' must be {ParticleType} or one of ('neutron', 'photon')."
+            )
 
         # convert particle type if needed
-        if isinstance(particle, str):
-            particle = ParticleType.from_string(particle)
+        particle = ParticleType(particle)
 
         if particle not in (ParticleType.NEUTRON, ParticleType.PHOTON):
             raise ValueError('Weight windows can only be applied for neutrons or photons')
diff --git a/openmc/particle_restart.py b/openmc/particle_restart.py
index 317d8761b61..44864107d9a 100644
--- a/openmc/particle_restart.py
+++ b/openmc/particle_restart.py
@@ -1,6 +1,7 @@
 import h5py
 
 import openmc.checkvalue as cv
+from .particle_type import ParticleType
 
 _VERSION_PARTICLE_RESTART = 2
 
@@ -28,8 +29,8 @@ class Particle:
         Type of simulation (criticality or fixed source)
     id : long
         Identifier of the particle
-    type : int
-        Particle type (1 = neutron, 2 = photon, 3 = electron, 4 = positron)
+    type : openmc.ParticleType
+        Particle type
     weight : float
         Weight of the particle
     energy : float
@@ -52,7 +53,7 @@ def __init__(self, filename):
             self.energy = f['energy'][()]
             self.generations_per_batch = f['generations_per_batch'][()]
             self.id = f['id'][()]
-            self.type = f['type'][()]
+            self.type = ParticleType(f['type'][()])
             self.n_particles = f['n_particles'][()]
             self.run_mode = f['run_mode'][()].decode()
             self.uvw = f['uvw'][()]
diff --git a/openmc/particle_type.py b/openmc/particle_type.py
new file mode 100644
index 00000000000..9d77134430a
--- /dev/null
+++ b/openmc/particle_type.py
@@ -0,0 +1,228 @@
+from numbers import Integral
+
+from openmc.data import gnds_name, zam, ATOMIC_SYMBOL
+
+
+_PDG_NAME = {
+    2112: 'neutron',
+    22: 'photon',
+    11: 'electron',
+    -11: 'positron',
+    2212: 'H1',
+}
+
+_ALIAS_PDG = {
+    'neutron': 2112,
+    'n': 2112,
+    'photon': 22,
+    'gamma': 22,
+    'electron': 11,
+    'positron': -11,
+    'proton': 2212,
+    'p': 2212,
+    'h1': 2212,
+    'deuteron': 1000010020,
+    'd': 1000010020,
+    'h2': 1000010020,
+    'triton': 1000010030,
+    't': 1000010030,
+    'h3': 1000010030,
+    'alpha': 1000020040,
+    'he4': 1000020040,
+}
+
+_LEGACY_PARTICLE_INDEX = {
+    0: 2112,
+    1: 22,
+    2: 11,
+    3: -11,
+}
+
+
+class ParticleType:
+    """Particle type defined by a PDG number.
+
+    ParticleType uses the Particle Data Group (PDG) Monte Carlo numbering scheme
+    to uniquely identify particle types. This includes elementary particles
+    (neutrons, photons, etc.) and nuclear codes for isotopes.
+
+    Parameters
+    ----------
+    value : str, int, or ParticleType
+        The particle identifier. Can be:
+
+        - A string name (e.g., 'neutron', 'photon', 'He4', 'U235')
+        - An integer PDG number (e.g., 2112 for neutron)
+        - A string with PDG prefix (e.g., 'pdg:2112')
+        - An existing ParticleType instance
+
+    Attributes
+    ----------
+    pdg_number : int
+        The PDG number for this particle type
+    zam : tuple of int or None
+        For nuclear particles, the (Z, A, m) tuple where Z is atomic number,
+        A is mass number, and m is metastable state. None for elementary particles.
+    is_nucleus : bool
+        Whether this particle is a nucleus (ion)
+
+    Examples
+    --------
+    >>> neutron = ParticleType('neutron')
+    >>> neutron.pdg_number
+    2112
+    >>> he4 = ParticleType('He4')
+    >>> he4.zam
+    (2, 4, 0)
+    >>> ParticleType(2112) == ParticleType('neutron')
+    True
+
+    """
+
+    __slots__ = ('_pdg_number',)
+
+    def __init__(self, value: 'str | int | ParticleType'):
+        if isinstance(value, ParticleType):
+            pdg = value._pdg_number
+        elif isinstance(value, str):
+            pdg = self._pdg_number_from_string(value)
+        elif isinstance(value, Integral):
+            pdg = int(value)
+            # Handle legacy particle indices (0, 1, 2, 3)
+            if pdg in _LEGACY_PARTICLE_INDEX:
+                pdg = _LEGACY_PARTICLE_INDEX[pdg]
+        else:
+            raise TypeError(f"Cannot create ParticleType from {type(value).__name__}")
+
+        self._pdg_number = pdg
+
+    def __eq__(self, other):
+        if isinstance(other, ParticleType):
+            return self._pdg_number == other._pdg_number
+        if isinstance(other, Integral):
+            return self._pdg_number == int(other)
+        if isinstance(other, str):
+            try:
+                return self._pdg_number == ParticleType(other)._pdg_number
+            except (ValueError, TypeError):
+                return False
+        return NotImplemented
+
+    def __hash__(self) -> int:
+        return hash(self._pdg_number)
+
+    def __int__(self) -> int:
+        return self._pdg_number
+
+    @property
+    def pdg_number(self) -> int:
+        return self._pdg_number
+
+    @staticmethod
+    def _pdg_number_from_string(value: str) -> int:
+        """Parse a string to get a PDG number.
+
+        Parameters
+        ----------
+        value : str
+            Particle identifier string
+
+        Returns
+        -------
+        int
+            PDG number
+
+        Raises
+        ------
+        ValueError
+            If string cannot be parsed as a valid particle identifier
+
+        """
+        s = value.strip()
+        if not s:
+            raise ValueError('Particle identifier cannot be empty.')
+
+        lower = s.lower()
+        if lower.startswith('pdg:'):
+            code_str = lower[4:]
+            try:
+                return int(code_str)
+            except ValueError:
+                raise ValueError(f'Invalid PDG number: {code_str}')
+
+        if lower in _ALIAS_PDG:
+            return _ALIAS_PDG[lower]
+
+        # Assume it is a GNDS nuclide name
+        Z, A, m = zam(s)
+        if Z <= 0 or Z > 999 or A <= 0 or A > 999 or m < 0 or m > 9:
+            raise ValueError('Invalid Z/A/m for nuclear PDG number.')
+        return 1000000000 + Z * 10000 + A * 10 + m
+
+    def __repr__(self) -> str:
+        return f'<ParticleType: {str(self)} (PDG={self._pdg_number})>'
+
+    def __str__(self) -> str:
+        """Return a canonical string representation of the particle type.
+
+        Returns
+        -------
+        str
+            Canonical name (e.g., 'neutron', 'He4', 'pdg:12345')
+
+        """
+        if self._pdg_number in _PDG_NAME:
+            return _PDG_NAME[self._pdg_number]
+
+        if (zam_tuple := self.zam) is not None:
+            Z, A, m = zam_tuple
+            if Z <= 0 or Z > max(ATOMIC_SYMBOL) or A <= 0 or A > 999:
+                raise ValueError(f"Invalid nuclear PDG number: {self._pdg_number}")
+            return gnds_name(Z, A, m)
+
+        return f'pdg:{self._pdg_number}'
+
+    @property
+    def zam(self) -> 'tuple[int, int, int] | None':
+        """Return the (Z, A, m) tuple for nuclear particles.
+
+        Returns
+        -------
+        tuple of int or None
+            For nuclear particles, returns (Z, A, m) where Z is atomic number,
+            A is mass number, and m is metastable state. Returns None for
+            elementary particles.
+
+        """
+        if self._pdg_number < 1000000000:
+            return None
+        Z = (self._pdg_number // 10000) % 1000
+        A = (self._pdg_number // 10) % 1000
+        m = self._pdg_number % 10
+        if Z <= 0 or A <= 0:
+            return None
+        else:
+            return (Z, A, m)
+
+    @property
+    def is_nucleus(self) -> bool:
+        """Return whether this particle is a nucleus.
+
+        Returns
+        -------
+        bool
+            True if the particle is a nucleus (ion), False otherwise
+
+        """
+        return self.zam is not None
+
+
+# Define common particle constants
+ParticleType.NEUTRON = ParticleType(2112)
+ParticleType.PHOTON = ParticleType(22)
+ParticleType.ELECTRON = ParticleType(11)
+ParticleType.POSITRON = ParticleType(-11)
+ParticleType.PROTON = ParticleType(2212)
+ParticleType.DEUTERON = ParticleType(1000010020)
+ParticleType.TRITON = ParticleType(1000010030)
+ParticleType.ALPHA = ParticleType(1000020040)
diff --git a/openmc/source.py b/openmc/source.py
index 5acf55c2194..a11cfd6c8d8 100644
--- a/openmc/source.py
+++ b/openmc/source.py
@@ -1,7 +1,6 @@
 from __future__ import annotations
 from abc import ABC, abstractmethod
 from collections.abc import Iterable, Sequence
-from enum import IntEnum
 from numbers import Real
 from pathlib import Path
 import warnings
@@ -19,6 +18,8 @@
 from openmc.stats.univariate import Univariate
 from ._xml import get_elem_list, get_text
 from .mesh import MeshBase, StructuredMesh, UnstructuredMesh
+from .particle_type import ParticleType
+from .statepoint import _VERSION_STATEPOINT
 from .utility_funcs import input_path
 
 
@@ -265,8 +266,8 @@ class IndependentSource(SourceBase):
         time distribution of source sites
     strength : float
         Strength of the source
-    particle : {'neutron', 'photon', 'electron', 'positron'}
-        Source particle type
+    particle : str or int or openmc.ParticleType
+        Source particle type (name, PDG number, or type)
     domains : iterable of openmc.Cell, openmc.Material, or openmc.Universe
         Domains to reject based on, i.e., if a sampled spatial location is not
         within one of these domains, it will be rejected.
@@ -302,10 +303,9 @@ class IndependentSource(SourceBase):
     type : str
         Indicator of source type: 'independent'
 
-    .. versionadded:: 0.14.0
-
-    particle : {'neutron', 'photon', 'electron', 'positron'}
-        Source particle type
+        .. versionadded:: 0.14.0
+    particle : str or int or openmc.ParticleType
+        Source particle type (alias, PDG number, or GNDS nuclide name)
     constraints : dict
         Constraints on sampled source particles. Valid keys include
         'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds',
@@ -320,7 +320,7 @@ def __init__(
         energy: openmc.stats.Univariate | None = None,
         time: openmc.stats.Univariate | None = None,
         strength: float = 1.0,
-        particle: str = 'neutron',
+        particle: str | int | ParticleType = 'neutron',
         domains: Sequence[openmc.Cell | openmc.Material |
                           openmc.Universe] | None = None,
         constraints: dict[str, Any] | None = None
@@ -405,14 +405,12 @@ def time(self, time):
         self._time = time
 
     @property
-    def particle(self):
+    def particle(self) -> ParticleType:
         return self._particle
 
     @particle.setter
     def particle(self, particle):
-        cv.check_value('source particle', particle,
-                       ['neutron', 'photon', 'electron', 'positron'])
-        self._particle = particle
+        self._particle = ParticleType(particle)
 
     def populate_xml_element(self, element):
         """Add necessary source information to an XML element
@@ -423,7 +421,7 @@ def populate_xml_element(self, element):
             XML element containing source data
 
         """
-        element.set("particle", self.particle)
+        element.set("particle", str(self.particle))
         if self.space is not None:
             element.append(self.space.to_xml_element())
         if self.angle is not None:
@@ -898,76 +896,6 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource:
         return cls(**kwargs)
 
 
-class ParticleType(IntEnum):
-    """
-    IntEnum class representing a particle type. Type
-    values mirror those found in the C++ class.
-    """
-    NEUTRON = 0
-    PHOTON = 1
-    ELECTRON = 2
-    POSITRON = 3
-
-    @classmethod
-    def from_string(cls, value: str):
-        """
-        Constructs a ParticleType instance from a string.
-
-        Parameters
-        ----------
-        value : str
-            The string representation of the particle type.
-
-        Returns
-        -------
-        The corresponding ParticleType instance.
-        """
-        try:
-            return cls[value.upper()]
-        except KeyError:
-            raise ValueError(
-                f"Invalid string for creation of {cls.__name__}: {value}")
-
-    @classmethod
-    def from_pdg_number(cls, pdg_number: int) -> ParticleType:
-        """Constructs a ParticleType instance from a PDG number.
-
-        The Particle Data Group at LBNL publishes a Monte Carlo particle
-        numbering scheme as part of the `Review of Particle Physics
-        <10.1103/PhysRevD.110.030001>`_. This method maps PDG numbers to the
-        corresponding :class:`ParticleType`.
-
-        Parameters
-        ----------
-        pdg_number : int
-            The PDG number of the particle type.
-
-        Returns
-        -------
-        The corresponding ParticleType instance.
-        """
-        try:
-            return {
-                2112: ParticleType.NEUTRON,
-                22: ParticleType.PHOTON,
-                11: ParticleType.ELECTRON,
-                -11: ParticleType.POSITRON,
-            }[pdg_number]
-        except KeyError:
-            raise ValueError(f"Unrecognized PDG number: {pdg_number}")
-
-    def __repr__(self) -> str:
-        """
-        Returns a string representation of the ParticleType instance.
-
-        Returns:
-            str: The lowercase name of the ParticleType instance.
-        """
-        return self.name.lower()
-
-    # needed for < Python 3.11
-    def __str__(self) -> str:
-        return self.__repr__()
 
 
 class SourceParticle:
@@ -992,8 +920,8 @@ class SourceParticle:
         Delayed group particle was created in (neutrons only)
     surf_id : int
         Surface ID where particle is at, if any.
-    particle : ParticleType
-        Type of the particle
+    particle : ParticleType or str or int
+        Type of the particle (type, name, or PDG number)
 
     """
 
@@ -1006,7 +934,7 @@ def __init__(
         wgt: float = 1.0,
         delayed_group: int = 0,
         surf_id: int = 0,
-        particle: ParticleType = ParticleType.NEUTRON
+        particle: ParticleType | str | int = ParticleType.NEUTRON
     ):
 
         self.r = tuple(r)
@@ -1018,9 +946,16 @@ def __init__(
         self.surf_id = surf_id
         self.particle = particle
 
+    @property
+    def particle(self) -> ParticleType:
+        return self._particle
+
+    @particle.setter
+    def particle(self, particle):
+        self._particle = ParticleType(particle)
+
     def __repr__(self):
-        name = self.particle.name.lower()
-        return f'<SourceParticle: {name} at E={self.E:.6e} eV>'
+        return f'<SourceParticle: {str(self.particle)} at E={self.E:.6e} eV>'
 
     def to_tuple(self) -> tuple:
         """Return source particle attributes as a tuple
@@ -1032,7 +967,7 @@ def to_tuple(self) -> tuple:
 
         """
         return (self.r, self.u, self.E, self.time, self.wgt,
-                self.delayed_group, self.surf_id, self.particle.value)
+                self.delayed_group, self.surf_id, self.particle.pdg_number)
 
 
 def write_source_file(
@@ -1116,12 +1051,7 @@ def from_mcpl(cls, filename: PathLike) -> ParticleList:
         particles = []
         with mcpl.MCPLFile(filename) as f:
             for particle in f.particles:
-                # Determine particle type based on the PDG number
-                try:
-                    particle_type = ParticleType.from_pdg_number(
-                        particle.pdgcode)
-                except ValueError:
-                    particle_type = "UNKNOWN"
+                particle_type = ParticleType(particle.pdgcode)
 
                 # Create a source particle instance. Note that MCPL stores
                 # energy in MeV and time in ms.
@@ -1179,7 +1109,7 @@ def to_dataframe(self) -> pd.DataFrame:
         # Extract the attributes of the source particles into a list of tuples
         data = [(sp.r[0], sp.r[1], sp.r[2], sp.u[0], sp.u[1], sp.u[2],
                  sp.E, sp.time, sp.wgt, sp.delayed_group, sp.surf_id,
-                 sp.particle.name.lower()) for sp in self]
+                 str(sp.particle)) for sp in self]
 
         # Define the column names for the DataFrame
         columns = ['x', 'y', 'z', 'u_x', 'u_y', 'u_z', 'E', 'time', 'wgt',
@@ -1226,6 +1156,7 @@ def export_to_hdf5(self, filename: PathLike, **kwargs):
         kwargs.setdefault('mode', 'w')
         with h5py.File(filename, **kwargs) as fh:
             fh.attrs['filetype'] = np.bytes_("source")
+            fh.attrs['version'] = np.array([_VERSION_STATEPOINT, 2])
             fh.create_dataset('source_bank', data=arr, dtype=source_dtype)
 
 
@@ -1337,7 +1268,7 @@ def read_collision_track_mcpl(file_path):
             data['material_id'].append(int(values_dict.get('material_id', 0)))
             data['universe_id'].append(int(values_dict.get('universe_id', 0)))
             data['n_collision'].append(int(values_dict.get('n_collision', 0)))
-            data['particle'].append(ParticleType.from_pdg_number(p.pdgcode))
+            data['particle'].append(ParticleType(p.pdgcode))
             data['parent_id'].append(int(values_dict.get('parent_id', 0)))
             data['progeny_id'].append(int(values_dict.get('progeny_id', 0)))
 
diff --git a/openmc/tracks.py b/openmc/tracks.py
index 81646e7d2e7..44e3c0eba6e 100644
--- a/openmc/tracks.py
+++ b/openmc/tracks.py
@@ -4,7 +4,8 @@
 import h5py
 
 from .checkvalue import check_filetype_version
-from .source import SourceParticle, ParticleType
+from .particle_type import ParticleType
+from .source import SourceParticle
 
 from pathlib import Path
 
@@ -25,7 +26,7 @@
 
 """
 def _particle_track_repr(self):
-    return f"<ParticleTrack: {self.particle}, {len(self.states)} states>"
+    return f"<ParticleTrack: {str(self.particle)}, {len(self.states)} states>"
 ParticleTrack.__repr__ = _particle_track_repr
 
 
@@ -92,8 +93,8 @@ def filter(self, particle=None, state_filter=None):
 
         Parameters
         ----------
-        particle : {'neutron', 'photon', 'electron', 'positron'}
-            Matching particle type
+        particle : str or int or openmc.ParticleType
+            Matching particle type (name, PDG number, or type)
         state_filter : function
             Function that takes a state (structured datatype) and returns a bool
             depending on some criteria.
@@ -126,7 +127,7 @@ def filter(self, particle=None, state_filter=None):
         for t in self:
             # Check for matching particle
             if particle is not None:
-                if t.particle.name.lower() != particle:
+                if t.particle != ParticleType(particle):
                     continue
 
             # Apply arbitrary state filter
@@ -184,7 +185,7 @@ def plot(self, axes=None):
     def sources(self):
         sources = []
         for particle_track in self:
-            particle_type = ParticleType(particle_track.particle)
+            particle_type = particle_track.particle
             state = particle_track.states[0]
             sources.append(
                 SourceParticle(
diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py
index 5d52a579a11..7797986df03 100644
--- a/openmc/weight_windows.py
+++ b/openmc/weight_windows.py
@@ -10,13 +10,12 @@
 import h5py
 
 import openmc
-from openmc.filter import _PARTICLES
 from openmc.mesh import MeshBase, RectilinearMesh, CylindricalMesh, SphericalMesh, UnstructuredMesh
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
 from ._xml import get_elem_list, get_text, clean_indentation
 from .mixin import IDManagerMixin
-from .utility_funcs import change_directory
+from .particle_type import ParticleType
 
 
 class WeightWindows(IDManagerMixin):
@@ -51,7 +50,7 @@ class WeightWindows(IDManagerMixin):
         A list of values for which each successive pair constitutes a range of
         energies in [eV] for a single bin. If no energy bins are provided, the
         maximum and minimum energy for the data available at runtime.
-    particle_type : {'neutron', 'photon'}
+    particle_type : str or int or openmc.ParticleType
         Particle type the weight windows apply to
     survival_ratio : float
         Ratio of the survival weight to the lower weight window bound for
@@ -116,7 +115,7 @@ def __init__(
         upper_ww_bounds: Iterable[float] | None = None,
         upper_bound_ratio: float | None = None,
         energy_bounds: Iterable[Real] | None = None,
-        particle_type: str = 'neutron',
+        particle_type: str | int | openmc.ParticleType = 'neutron',
         survival_ratio: float = 3.0,
         max_lower_bound_ratio: float | None = None,
         max_split: int = 10,
@@ -213,13 +212,15 @@ def mesh(self, mesh: MeshBase):
         self._mesh = mesh
 
     @property
-    def particle_type(self) -> str:
+    def particle_type(self) -> ParticleType:
         return self._particle_type
 
     @particle_type.setter
-    def particle_type(self, pt: str):
-        cv.check_value('Particle type', pt, _PARTICLES)
-        self._particle_type = pt
+    def particle_type(self, pt):
+        ptype = ParticleType(pt)
+        if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}:
+            raise ValueError("Weight windows can only be applied for neutrons or photons")
+        self._particle_type = ptype
 
     @property
     def energy_bounds(self) -> Iterable[Real]:
@@ -329,7 +330,7 @@ def to_xml_element(self) -> ET.Element:
         subelement.text = str(self.mesh.id)
 
         subelement = ET.SubElement(element, 'particle_type')
-        subelement.text = self.particle_type
+        subelement.text = str(self.particle_type)
 
         if self.energy_bounds is not None:
             subelement = ET.SubElement(element, 'energy_bounds')
@@ -494,7 +495,7 @@ class WeightWindowGenerator:
         A list of values for which each successive pair constitutes a range of
         energies in [eV] for a single bin. If no energy bins are provided, the
         maximum and minimum energy for the data available at runtime.
-    particle_type : {'neutron', 'photon'}
+    particle_type : str or int or openmc.ParticleType
         Particle type the weight windows apply to
     method : {'magic', 'fw_cadis'}
         The weight window generation methodology applied during an update.
@@ -513,7 +514,7 @@ class WeightWindowGenerator:
     energy_bounds : Iterable of Real
         A list of values for which each successive pair constitutes a range of
         energies in [eV] for a single bin
-    particle_type : {'neutron', 'photon'}
+    particle_type : openmc.ParticleType
         Particle type the weight windows apply to
     method : {'magic', 'fw_cadis'}
         The weight window generation methodology applied during an update.
@@ -534,7 +535,7 @@ def __init__(
         self,
         mesh: openmc.MeshBase,
         energy_bounds: Sequence[float] | None = None,
-        particle_type: str = 'neutron',
+        particle_type: str | int | openmc.ParticleType = 'neutron',
         method: str = 'magic',
         max_realizations: int = 1,
         update_interval: int = 1,
@@ -555,7 +556,7 @@ def __init__(
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += f'\t{"Mesh":<20}=\t{self.mesh.id}\n'
-        string += f'\t{"Particle:":<20}=\t{self.particle_type}\n'
+        string += f'\t{"Particle:":<20}=\t{str(self.particle_type)}\n'
         string += f'\t{"Energy Bounds:":<20}=\t{self.energy_bounds}\n'
         string += f'\t{"Method":<20}=\t{self.method}\n'
         string += f'\t{"Max Realizations:":<20}=\t{self.max_realizations}\n'
@@ -586,13 +587,15 @@ def energy_bounds(self, eb: Iterable[float]):
         self._energy_bounds = eb
 
     @property
-    def particle_type(self) -> str:
+    def particle_type(self) -> ParticleType:
         return self._particle_type
 
     @particle_type.setter
-    def particle_type(self, pt: str):
-        cv.check_value('particle type', pt, ('neutron', 'photon'))
-        self._particle_type = pt
+    def particle_type(self, pt):
+        ptype = ParticleType(pt)
+        if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}:
+            raise ValueError("Weight windows can only be applied for neutrons or photons")
+        self._particle_type = ptype
 
     @property
     def method(self) -> str:
@@ -695,7 +698,7 @@ def to_xml_element(self):
             subelement = ET.SubElement(element, 'energy_bounds')
             subelement.text = ' '.join(str(e) for e in self.energy_bounds)
         particle_elem = ET.SubElement(element, 'particle_type')
-        particle_elem.text = self.particle_type
+        particle_elem.text = str(self.particle_type)
         realizations_elem = ET.SubElement(element, 'max_realizations')
         realizations_elem.text = str(self.max_realizations)
         update_interval_elem = ET.SubElement(element, 'update_interval')
@@ -730,7 +733,7 @@ def from_xml_element(cls, elem: ET.Element, meshes: dict) -> Self:
 
         mesh_id = int(get_text(elem, 'mesh'))
         mesh = meshes[mesh_id]
-        
+
         energy_bounds = get_elem_list(elem, "energy_bounds, float")
         particle_type = get_text(elem, 'particle_type')
 
diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp
index a2320e0b469..d77066fb0eb 100644
--- a/src/bremsstrahlung.cpp
+++ b/src/bremsstrahlung.cpp
@@ -31,13 +31,13 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
   if (p.material() == MATERIAL_VOID)
     return;
 
-  int photon = static_cast<int>(ParticleType::photon);
+  int photon = ParticleType::photon().transport_index();
   if (p.E() < settings::energy_cutoff[photon])
     return;
 
   // Get bremsstrahlung data for this material and particle type
   BremsstrahlungData* mat;
-  if (p.type() == ParticleType::positron) {
+  if (p.type() == ParticleType::positron()) {
     mat = &model::materials[p.material()]->ttb_->positron;
   } else {
     mat = &model::materials[p.material()]->ttb_->electron;
@@ -119,7 +119,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
       }
 
       // Create secondary photon
-      p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon);
+      p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon());
       *E_lost += w;
     }
   }
diff --git a/src/material.cpp b/src/material.cpp
index 54caa38409a..072e6decad1 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -819,9 +819,9 @@ void Material::calculate_xs(Particle& p) const
   p.macro_xs().fission = 0.0;
   p.macro_xs().nu_fission = 0.0;
 
-  if (p.type() == ParticleType::neutron) {
+  if (p.type().is_neutron()) {
     this->calculate_neutron_xs(p);
-  } else if (p.type() == ParticleType::photon) {
+  } else if (p.type().is_photon()) {
     this->calculate_photon_xs(p);
   }
 }
@@ -829,7 +829,7 @@ void Material::calculate_xs(Particle& p) const
 void Material::calculate_neutron_xs(Particle& p) const
 {
   // Find energy index on energy grid
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   int i_grid =
     std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing;
 
diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp
index 256f3343fc6..30c41ec5afa 100644
--- a/src/mcpl_interface.cpp
+++ b/src/mcpl_interface.cpp
@@ -321,25 +321,7 @@ inline void ensure_mcpl_ready_or_fatal()
 SourceSite mcpl_particle_to_site(const mcpl_particle_repr_t* particle_repr)
 {
   SourceSite site;
-  switch (particle_repr->pdgcode) {
-  case 2112:
-    site.particle = ParticleType::neutron;
-    break;
-  case 22:
-    site.particle = ParticleType::photon;
-    break;
-  case 11:
-    site.particle = ParticleType::electron;
-    break;
-  case -11:
-    site.particle = ParticleType::positron;
-    break;
-  default:
-    fatal_error(fmt::format(
-      "MCPL: Encountered unexpected PDG code {} when converting to SourceSite.",
-      particle_repr->pdgcode));
-    break;
-  }
+  site.particle = ParticleType {particle_repr->pdgcode};
 
   // Copy position and direction
   site.r.x = particle_repr->position[0];
@@ -368,7 +350,6 @@ vector<SourceSite> mcpl_source_sites(std::string path)
   }
 
   size_t n_particles_in_file = g_mcpl_api->hdr_nparticles(mcpl_file);
-  size_t n_skipped = 0;
   if (n_particles_in_file > 0) {
     sites.reserve(n_particles_in_file);
   }
@@ -381,31 +362,16 @@ vector<SourceSite> mcpl_source_sites(std::string path)
         path, sites.size(), n_particles_in_file));
       break;
     }
-    if (p_repr->pdgcode == 2112 || p_repr->pdgcode == 22 ||
-        p_repr->pdgcode == 11 || p_repr->pdgcode == -11) {
-      sites.push_back(mcpl_particle_to_site(p_repr));
-    } else {
-      n_skipped++;
-    }
+    sites.push_back(mcpl_particle_to_site(p_repr));
   }
 
   g_mcpl_api->close_file(mcpl_file);
 
-  if (n_skipped > 0 && n_particles_in_file > 0) {
-    double percent_skipped =
-      100.0 * static_cast<double>(n_skipped) / n_particles_in_file;
-    warning(fmt::format(
-      "MCPL: Skipped {} of {} total particles ({:.1f}%) in file '{}' because "
-      "their type is not supported by OpenMC.",
-      n_skipped, n_particles_in_file, percent_skipped, path));
-  }
-
   if (sites.empty()) {
     if (n_particles_in_file > 0) {
       fatal_error(fmt::format(
-        "MCPL file '{}' contained {} particles, but none were of the supported "
-        "types (neutron, photon, electron, positron). OpenMC cannot proceed "
-        "without source particles.",
+        "MCPL file '{}' contained {} particles, but no particles could be "
+        "read.",
         path, n_particles_in_file));
     } else {
       fatal_error(fmt::format(
@@ -461,22 +427,7 @@ void write_mcpl_source_bank_internal(mcpl_outfile_t* file_id,
         p_repr.ekin = site.E * 1e-6;
         p_repr.time = site.time * 1e3;
         p_repr.weight = site.wgt;
-        switch (site.particle) {
-        case ParticleType::neutron:
-          p_repr.pdgcode = 2112;
-          break;
-        case ParticleType::photon:
-          p_repr.pdgcode = 22;
-          break;
-        case ParticleType::electron:
-          p_repr.pdgcode = 11;
-          break;
-        case ParticleType::positron:
-          p_repr.pdgcode = -11;
-          break;
-        default:
-          continue;
-        }
+        p_repr.pdgcode = site.particle.pdg_number();
         g_mcpl_api->add_particle(file_id, &p_repr);
       }
     }
@@ -633,22 +584,7 @@ void write_mcpl_collision_track_internal(mcpl_outfile_t* file_id,
       p_repr.ekin = site.E * 1e-6;
       p_repr.time = site.time * 1e3;
       p_repr.weight = site.wgt;
-      switch (site.particle) {
-      case ParticleType::neutron:
-        p_repr.pdgcode = 2112;
-        break;
-      case ParticleType::photon:
-        p_repr.pdgcode = 22;
-        break;
-      case ParticleType::electron:
-        p_repr.pdgcode = 11;
-        break;
-      case ParticleType::positron:
-        p_repr.pdgcode = -11;
-        break;
-      default:
-        continue;
-      }
+      p_repr.pdgcode = site.particle.pdg_number();
       g_mcpl_api->add_particle(file_id, &p_repr);
     }
   } else {
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 6ee15b9c18b..a48530caefa 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -376,7 +376,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
 
     SourceSite site;
     site.E = 1.0;
-    site.particle = ParticleType::neutron;
+    site.particle = ParticleType::neutron();
 
     for (int axis = 0; axis < 3; ++axis) {
       // Set starting position and direction
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index ed734d401ea..34f87d17984 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -244,7 +244,7 @@ void MgxsInterface::read_header(const std::string& path_cross_sections)
 void put_mgxs_header_data_to_globals()
 {
   // Get the minimum and maximum energies
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   data::energy_min[neutron] = data::mg.energy_bins_.back();
   data::energy_max[neutron] = data::mg.energy_bins_.front();
 
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 5ae6e30ee24..69e603a7c66 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -379,7 +379,7 @@ void Nuclide::create_derived(
       auto pprod = xt::view(xs_[t], xt::range(j, j + n), XS_PHOTON_PROD);
 
       for (const auto& p : rx->products_) {
-        if (p.particle_ == ParticleType::photon) {
+        if (p.particle_.is_photon()) {
           for (int k = 0; k < n; ++k) {
             double E = grid_[t].energy[k + j];
 
@@ -501,7 +501,7 @@ void Nuclide::create_derived(
 
 void Nuclide::init_grid()
 {
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   double E_min = data::energy_min[neutron];
   double E_max = data::energy_max[neutron];
   int M = settings::n_log_bins;
@@ -552,7 +552,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
         for (int i = 1; i < rx->products_.size(); ++i) {
           // Skip any non-neutron products
           const auto& product = rx->products_[i];
-          if (product.particle_ != ParticleType::neutron)
+          if (!product.particle_.is_neutron())
             continue;
 
           // Evaluate yield
diff --git a/src/output.cpp b/src/output.cpp
index 80e2b10ab89..ae2daaffc14 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -155,21 +155,21 @@ std::string time_stamp()
 void print_particle(Particle& p)
 {
   // Display particle type and ID.
-  switch (p.type()) {
-  case ParticleType::neutron:
+  switch (p.type().pdg_number()) {
+  case PDG_NEUTRON:
     fmt::print("Neutron ");
     break;
-  case ParticleType::photon:
+  case PDG_PHOTON:
     fmt::print("Photon ");
     break;
-  case ParticleType::electron:
+  case PDG_ELECTRON:
     fmt::print("Electron ");
     break;
-  case ParticleType::positron:
+  case PDG_POSITRON:
     fmt::print("Positron ");
     break;
   default:
-    fmt::print("Unknown Particle ");
+    fmt::print("Particle {} ", p.type().str());
   }
   fmt::print("{}\n", p.id());
 
diff --git a/src/particle.cpp b/src/particle.cpp
index b9f9c8f86b6..a1176abc79a 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -48,17 +48,19 @@ double Particle::speed() const
   if (settings::run_CE) {
     // Determine mass in eV/c^2
     double mass;
-    switch (this->type()) {
-    case ParticleType::neutron:
+    switch (this->type().pdg_number()) {
+    case PDG_NEUTRON:
       mass = MASS_NEUTRON_EV;
       break;
-    case ParticleType::photon:
+    case PDG_PHOTON:
       mass = 0.0;
       break;
-    case ParticleType::electron:
-    case ParticleType::positron:
+    case PDG_ELECTRON:
+    case PDG_POSITRON:
       mass = MASS_ELECTRON_EV;
       break;
+    default:
+      fatal_error("Unsupported particle for speed calculation.");
     }
     // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
     return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
@@ -77,7 +79,11 @@ bool Particle::create_secondary(
 {
   // If energy is below cutoff for this particle, don't create secondary
   // particle
-  if (E < settings::energy_cutoff[static_cast<int>(type)]) {
+  int idx = type.transport_index();
+  if (idx == C_NONE) {
+    return false;
+  }
+  if (E < settings::energy_cutoff[idx]) {
     return false;
   }
 
@@ -235,7 +241,8 @@ void Particle::event_advance()
   boundary() = distance_to_boundary(*this);
 
   // Sample a distance to collision
-  if (type() == ParticleType::electron || type() == ParticleType::positron) {
+  if (type() == ParticleType::electron() ||
+      type() == ParticleType::positron()) {
     collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
   } else if (macro_xs().total == 0.0) {
     collision_distance() = INFINITY;
@@ -244,7 +251,7 @@ void Particle::event_advance()
   }
 
   double speed = this->speed();
-  double time_cutoff = settings::time_cutoff[static_cast<int>(type())];
+  double time_cutoff = settings::time_cutoff[type().transport_index()];
   double distance_cutoff =
     (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
 
@@ -269,8 +276,7 @@ void Particle::event_advance()
   }
 
   // Score track-length estimate of k-eff
-  if (settings::run_mode == RunMode::EIGENVALUE &&
-      type() == ParticleType::neutron) {
+  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
     keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
   }
 
@@ -331,8 +337,7 @@ void Particle::event_cross_surface()
 void Particle::event_collide()
 {
   // Score collision estimate of keff
-  if (settings::run_mode == RunMode::EIGENVALUE &&
-      type() == ParticleType::neutron) {
+  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
     keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
   }
 
@@ -370,8 +375,7 @@ void Particle::event_collide()
     }
   }
 
-  if (!model::active_pulse_height_tallies.empty() &&
-      type() == ParticleType::photon) {
+  if (!model::active_pulse_height_tallies.empty() && type().is_photon()) {
     pht_collision_energy();
   }
 
@@ -442,7 +446,7 @@ void Particle::event_revive_from_secondary()
 
     // Subtract secondary particle energy from interim pulse-height results
     if (!model::active_pulse_height_tallies.empty() &&
-        this->type() == ParticleType::photon) {
+        this->type().is_photon()) {
       // Since the birth cell of the particle has not been set we
       // have to determine it before the energy of the secondary particle can be
       // removed from the pulse-height of this cell.
@@ -525,7 +529,7 @@ void Particle::pht_collision_energy()
 
     // If the energy of the particle is below the cutoff, it will not be sampled
     // so its energy is added to the pulse-height in the cell
-    int photon = static_cast<int>(ParticleType::photon);
+    int photon = ParticleType::photon().transport_index();
     if (E() < settings::energy_cutoff[photon]) {
       pht_storage()[index] += E();
     }
@@ -824,7 +828,7 @@ void Particle::write_restart() const
       break;
     }
     write_dataset(file_id, "id", id());
-    write_dataset(file_id, "type", static_cast<int>(type()));
+    write_dataset(file_id, "type", type().pdg_number());
 
     int64_t i = current_work();
     if (settings::run_mode == RunMode::EIGENVALUE) {
@@ -878,37 +882,6 @@ void Particle::update_neutron_xs(
 //==============================================================================
 // Non-method functions
 //==============================================================================
-
-std::string particle_type_to_str(ParticleType type)
-{
-  switch (type) {
-  case ParticleType::neutron:
-    return "neutron";
-  case ParticleType::photon:
-    return "photon";
-  case ParticleType::electron:
-    return "electron";
-  case ParticleType::positron:
-    return "positron";
-  }
-  UNREACHABLE();
-}
-
-ParticleType str_to_particle_type(std::string str)
-{
-  if (str == "neutron") {
-    return ParticleType::neutron;
-  } else if (str == "photon") {
-    return ParticleType::photon;
-  } else if (str == "electron") {
-    return ParticleType::electron;
-  } else if (str == "positron") {
-    return ParticleType::positron;
-  } else {
-    throw std::invalid_argument {fmt::format("Invalid particle name: {}", str)};
-  }
-}
-
 void add_surf_source_to_bank(Particle& p, const Surface& surf)
 {
   if (simulation::current_batch <= settings::n_inactive ||
diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp
index 6b7778211af..f02fcb94b55 100644
--- a/src/particle_restart.cpp
+++ b/src/particle_restart.cpp
@@ -31,6 +31,16 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   hid_t file_id = file_open(settings::path_particle_restart, 'r');
 
   // Read data from file
+  bool legacy_particle_codes = true;
+  if (attribute_exists(file_id, "version")) {
+    array<int, 2> version;
+    read_attribute(file_id, "version", version);
+    if (version[0] > VERSION_PARTICLE_RESTART[0] ||
+        (version[0] == VERSION_PARTICLE_RESTART[0] && version[1] >= 1)) {
+      legacy_particle_codes = false;
+    }
+  }
+
   read_dataset(file_id, "current_batch", simulation::current_batch);
   read_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
   read_dataset(file_id, "current_generation", simulation::current_gen);
@@ -45,7 +55,8 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
   read_dataset(file_id, "id", p.id());
   int type;
   read_dataset(file_id, "type", type);
-  p.type() = static_cast<ParticleType>(type);
+  p.type() = legacy_particle_codes ? legacy_particle_index_to_type(type)
+                                   : ParticleType {type};
   read_dataset(file_id, "weight", p.wgt());
   read_dataset(file_id, "energy", p.E());
   read_dataset(file_id, "xyz", p.r());
diff --git a/src/particle_type.cpp b/src/particle_type.cpp
new file mode 100644
index 00000000000..fe8f9fe2a01
--- /dev/null
+++ b/src/particle_type.cpp
@@ -0,0 +1,246 @@
+#include "openmc/particle_type.h"
+
+#include <algorithm>
+#include <cctype>
+#include <stdexcept>
+
+#include "openmc/string_utils.h"
+
+namespace openmc {
+namespace {
+
+constexpr const char* ATOMIC_SYMBOL[] = {"", "H", "He", "Li", "Be", "B", "C",
+  "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca",
+  "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As",
+  "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd",
+  "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr",
+  "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf",
+  "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At",
+  "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf",
+  "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg",
+  "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"};
+
+constexpr int MAX_Z =
+  static_cast<int>(sizeof(ATOMIC_SYMBOL) / sizeof(ATOMIC_SYMBOL[0])) - 1;
+
+bool is_integer_string(const std::string& s)
+{
+  if (s.empty())
+    return false;
+  size_t i = 0;
+  if (s[0] == '-' || s[0] == '+') {
+    if (s.size() == 1)
+      return false;
+    i = 1;
+  }
+  for (; i < s.size(); ++i) {
+    if (!std::isdigit(static_cast<unsigned char>(s[i])))
+      return false;
+  }
+  return true;
+}
+
+int atomic_number_from_symbol(std::string_view symbol)
+{
+  for (int z = 1; z <= MAX_Z; ++z) {
+    if (symbol == ATOMIC_SYMBOL[z]) {
+      return z;
+    }
+  }
+  return 0;
+}
+
+bool parse_gnds_nuclide(std::string_view name, int& Z, int& A, int& m)
+{
+  if (name.empty())
+    return false;
+
+  size_t pos = 0;
+  if (!std::isupper(static_cast<unsigned char>(name[pos])))
+    return false;
+
+  std::string symbol;
+  symbol += name[pos++];
+  if (pos < name.size() &&
+      std::islower(static_cast<unsigned char>(name[pos]))) {
+    symbol += name[pos++];
+  }
+
+  if (pos >= name.size() ||
+      !std::isdigit(static_cast<unsigned char>(name[pos]))) {
+    return false;
+  }
+
+  size_t a_start = pos;
+  while (
+    pos < name.size() && std::isdigit(static_cast<unsigned char>(name[pos]))) {
+    ++pos;
+  }
+  A = std::stoi(std::string {name.substr(a_start, pos - a_start)});
+  if (A <= 0 || A > 999)
+    return false;
+
+  m = 0;
+  if (pos < name.size()) {
+    if (name[pos] != '_' || pos + 2 >= name.size() || name[pos + 1] != 'm') {
+      return false;
+    }
+    pos += 2;
+    size_t m_start = pos;
+    while (pos < name.size() &&
+           std::isdigit(static_cast<unsigned char>(name[pos]))) {
+      ++pos;
+    }
+    if (m_start == pos)
+      return false;
+    m = std::stoi(std::string {name.substr(m_start, pos - m_start)});
+    if (m < 0 || m > 9)
+      return false;
+  }
+
+  if (pos != name.size())
+    return false;
+
+  Z = atomic_number_from_symbol(symbol);
+  return Z != 0;
+}
+
+// Helper to convert nuclear PDG number to nuclide name
+std::string nuclide_name_from_pdg(int32_t pdg)
+{
+  int32_t code = pdg;
+  int m = code % 10;
+  int A = (code / 10) % 1000;
+  int Z = (code / 10000) % 1000;
+
+  if (Z <= 0 || Z > MAX_Z || A <= 0 || A > 999) {
+    throw std::invalid_argument {
+      "Invalid nuclear PDG number: " + std::to_string(pdg)};
+  }
+
+  std::string name = ATOMIC_SYMBOL[Z] + std::to_string(A);
+  if (m > 0) {
+    name += "_m" + std::to_string(m);
+  }
+  return name;
+}
+
+} // namespace
+
+//==============================================================================
+// ParticleType member function implementations
+//==============================================================================
+
+ParticleType::ParticleType(std::string_view str)
+{
+  std::string s {str};
+  strtrim(s);
+  if (s.empty()) {
+    throw std::invalid_argument {"Particle string is empty."};
+  }
+
+  std::string lower = s;
+  to_lower(lower);
+
+  // Check for pdg: prefix
+  if (starts_with(lower, "pdg:")) {
+    std::string value_str = lower.substr(4);
+    if (!is_integer_string(value_str)) {
+      throw std::invalid_argument {"Invalid PDG number: " + value_str};
+    }
+    pdg_number_ = std::stoi(value_str);
+    return;
+  }
+
+  // Check for known particle names
+  if (lower == "neutron" || lower == "n") {
+    pdg_number_ = PDG_NEUTRON;
+    return;
+  }
+  if (lower == "photon" || lower == "gamma") {
+    pdg_number_ = PDG_PHOTON;
+    return;
+  }
+  if (lower == "electron") {
+    pdg_number_ = PDG_ELECTRON;
+    return;
+  }
+  if (lower == "positron") {
+    pdg_number_ = PDG_POSITRON;
+    return;
+  }
+  if (lower == "proton" || lower == "p" || lower == "h1") {
+    pdg_number_ = PDG_PROTON;
+    return;
+  }
+  if (lower == "deuteron" || lower == "d" || lower == "h2") {
+    pdg_number_ = PDG_DEUTERON;
+    return;
+  }
+  if (lower == "triton" || lower == "t" || lower == "h3") {
+    pdg_number_ = PDG_TRITON;
+    return;
+  }
+  if (lower == "alpha" || lower == "he4") {
+    pdg_number_ = PDG_ALPHA;
+    return;
+  }
+
+  // Check for integer string
+  if (is_integer_string(s)) {
+    pdg_number_ = std::stoi(s);
+    return;
+  }
+
+  // Try to parse as GNDS nuclide name
+  int Z = 0;
+  int A = 0;
+  int m = 0;
+  if (!parse_gnds_nuclide(s, Z, A, m)) {
+    throw std::invalid_argument {"Invalid nuclide name: " + s};
+  }
+  pdg_number_ = 1000000000 + Z * 10000 + A * 10 + m;
+}
+
+std::string ParticleType::str() const
+{
+  if (pdg_number_ == PDG_NEUTRON)
+    return "neutron";
+  if (pdg_number_ == PDG_PHOTON)
+    return "photon";
+  if (pdg_number_ == PDG_ELECTRON)
+    return "electron";
+  if (pdg_number_ == PDG_POSITRON)
+    return "positron";
+  if (pdg_number_ == PDG_PROTON)
+    return "proton";
+
+  if (is_nucleus()) {
+    return nuclide_name_from_pdg(pdg_number_);
+  }
+
+  return "pdg:" + std::to_string(pdg_number_);
+}
+
+//==============================================================================
+// Free function implementations
+//==============================================================================
+
+ParticleType legacy_particle_index_to_type(int index)
+{
+  switch (index) {
+  case 0:
+    return ParticleType {PDG_NEUTRON};
+  case 1:
+    return ParticleType {PDG_PHOTON};
+  case 2:
+    return ParticleType {PDG_ELECTRON};
+  case 3:
+    return ParticleType {PDG_POSITRON};
+  default:
+    throw std::invalid_argument {
+      "Invalid legacy particle index: " + std::to_string(index)};
+  }
+}
+
+} // namespace openmc
diff --git a/src/photon.cpp b/src/photon.cpp
index 4926e3eae6a..951acb9fbd8 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -293,7 +293,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
     close_group(rgroup);
 
     // Truncate the bremsstrahlung data at the cutoff energy
-    int photon = static_cast<int>(ParticleType::photon);
+    int photon = ParticleType::photon().transport_index();
     const auto& E {electron_energy};
     double cutoff = settings::energy_cutoff[photon];
     if (cutoff > E(0)) {
@@ -805,7 +805,7 @@ void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const
     if (shell.transitions.empty()) {
       Direction u = isotropic_direction(p.current_seed());
       double E = shell.binding_energy;
-      p.create_secondary(p.wgt(), u, E, ParticleType::photon);
+      p.create_secondary(p.wgt(), u, E, ParticleType::photon());
       continue;
     }
 
@@ -833,12 +833,13 @@ void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const
       holes[n_holes++] = transition.secondary_subshell;
 
       // Create auger electron
-      p.create_secondary(p.wgt(), u, transition.energy, ParticleType::electron);
+      p.create_secondary(
+        p.wgt(), u, transition.energy, ParticleType::electron());
     } else {
       // Radiative transition -- get X-ray energy
 
       // Create fluorescent photon
-      p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon);
+      p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon());
     }
   }
 }
diff --git a/src/physics.cpp b/src/physics.cpp
index 41509af97be..05a9e59d940 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -46,27 +46,29 @@ void collision(Particle& p)
   ++(p.n_collision());
 
   // Sample reaction for the material the particle is in
-  switch (p.type()) {
-  case ParticleType::neutron:
+  switch (p.type().pdg_number()) {
+  case PDG_NEUTRON:
     sample_neutron_reaction(p);
     break;
-  case ParticleType::photon:
+  case PDG_PHOTON:
     sample_photon_reaction(p);
     break;
-  case ParticleType::electron:
+  case PDG_ELECTRON:
     sample_electron_reaction(p);
     break;
-  case ParticleType::positron:
+  case PDG_POSITRON:
     sample_positron_reaction(p);
     break;
+  default:
+    fatal_error("Unsupported particle PDG for collision sampling.");
   }
 
   if (settings::weight_window_checkpoint_collision)
     apply_weight_windows(p);
 
   // Kill particle if energy falls below cutoff
-  int type = static_cast<int>(p.type());
-  if (p.E() < settings::energy_cutoff[type]) {
+  int type = p.type().transport_index();
+  if (type != C_NONE && p.E() < settings::energy_cutoff[type]) {
     p.wgt() = 0.0;
   }
 
@@ -75,11 +77,11 @@ void collision(Particle& p)
     std::string msg;
     if (p.event() == TallyEvent::KILL) {
       msg = fmt::format("    Killed. Energy = {} eV.", p.E());
-    } else if (p.type() == ParticleType::neutron) {
+    } else if (p.type().is_neutron()) {
       msg = fmt::format("    {} with {}. Energy = {} eV.",
         reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_,
         p.E());
-    } else if (p.type() == ParticleType::photon) {
+    } else if (p.type().is_photon()) {
       msg = fmt::format("    {} with {}. Energy = {} eV.",
         reaction_name(p.event_mt()),
         to_element(data::nuclides[p.event_nuclide()]->name_), p.E());
@@ -208,7 +210,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     // Initialize fission site object with particle data
     SourceSite site;
     site.r = p.r();
-    site.particle = ParticleType::neutron;
+    site.particle = ParticleType::neutron();
     site.time = p.time();
     site.wgt = 1. / weight;
     site.surf_id = 0;
@@ -218,7 +220,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
 
     // Reject site if it exceeds time cutoff
     if (site.delayed_group > 0) {
-      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
       if (site.time > t_cutoff) {
         continue;
       }
@@ -289,7 +291,7 @@ void sample_photon_reaction(Particle& p)
   // Kill photon if below energy cutoff -- an extra check is made here because
   // photons with energy below the cutoff may have been produced by neutrons
   // reactions or atomic relaxation
-  int photon = static_cast<int>(ParticleType::photon);
+  int photon = ParticleType::photon().transport_index();
   if (p.E() < settings::energy_cutoff[photon]) {
     p.E() = 0.0;
     p.wgt() = 0.0;
@@ -339,13 +341,13 @@ void sample_photon_reaction(Particle& p)
     // Create Compton electron
     double phi = uniform_distribution(0., 2.0 * PI, p.current_seed());
     double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b;
-    int electron = static_cast<int>(ParticleType::electron);
+    int electron = ParticleType::electron().transport_index();
     if (E_electron >= settings::energy_cutoff[electron]) {
       double mu_electron = (alpha - alpha_out * p.mu()) /
                            std::sqrt(alpha * alpha + alpha_out * alpha_out -
                                      2.0 * alpha * alpha_out * p.mu());
       Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed());
-      p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
+      p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
     }
 
     // Allow electrons to fill orbital and produce Auger electrons and
@@ -418,7 +420,7 @@ void sample_photon_reaction(Particle& p)
         u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi);
 
         // Create secondary electron
-        p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
+        p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
 
         // Allow electrons to fill orbital and produce auger electrons
         // and fluorescent photons
@@ -443,12 +445,11 @@ void sample_photon_reaction(Particle& p)
 
     // Create secondary electron
     Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed());
-    p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
+    p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
 
     // Create secondary positron
     u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed());
-    p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron);
-
+    p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron());
     p.event() = TallyEvent::ABSORB;
     p.event_mt() = PAIR_PROD;
     p.wgt() = 0.0;
@@ -483,8 +484,8 @@ void sample_positron_reaction(Particle& p)
   Direction u = isotropic_direction(p.current_seed());
 
   // Create annihilation photon pair traveling in opposite directions
-  p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon);
-  p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon);
+  p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon());
+  p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon());
 
   p.E() = 0.0;
   p.wgt() = 0.0;
@@ -609,7 +610,7 @@ void sample_photon_product(
       continue;
 
     for (int j = 0; j < rx->products_.size(); ++j) {
-      if (rx->products_[j].particle_ == ParticleType::photon) {
+      if (rx->products_[j].particle_.is_photon()) {
         // For fission, artificially increase the photon yield to account
         // for delayed photons
         double f = 1.0;
@@ -1098,7 +1099,7 @@ void sample_fission_neutron(
     rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed);
 
     // resample if energy is greater than maximum neutron energy
-    constexpr int neutron = static_cast<int>(ParticleType::neutron);
+    int neutron = ParticleType::neutron().transport_index();
     if (site->E < data::energy_max[neutron])
       break;
 
@@ -1158,7 +1159,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   if (std::floor(yield) == yield && yield > 0) {
     // If yield is integral, create exactly that many secondary particles
     for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
-      p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron);
+      p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron());
     }
   } else {
     // Otherwise, change weight of particle based on yield
@@ -1214,7 +1215,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
     }
 
     // Create the secondary photon
-    bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon);
+    bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon());
 
     // Tag secondary particle with parent nuclide
     if (created_photon && settings::use_decay_photons) {
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 4c28cb1795a..866d4d728e8 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -136,7 +136,7 @@ void create_fission_sites(Particle& p)
     // Initialize fission site object with particle data
     SourceSite site;
     site.r = p.r();
-    site.particle = ParticleType::neutron;
+    site.particle = ParticleType::neutron();
     site.time = p.time();
     site.wgt = 1. / weight;
 
@@ -171,7 +171,7 @@ void create_fission_sites(Particle& p)
       site.time -= std::log(prn(p.current_seed())) / decay_rate;
 
       // Reject site if it exceeds time cutoff
-      double t_cutoff = settings::time_cutoff[static_cast<int>(site.particle)];
+      double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
       if (site.time > t_cutoff) {
         continue;
       }
diff --git a/src/reaction.cpp b/src/reaction.cpp
index d96790c6d43..9ac5a1f528a 100644
--- a/src/reaction.cpp
+++ b/src/reaction.cpp
@@ -70,9 +70,8 @@ Reaction::Reaction(
 
   if (settings::use_decay_photons) {
     // Remove photon products for D1S method
-    products_.erase(
-      std::remove_if(products_.begin(), products_.end(),
-        [](const auto& p) { return p.particle_ == ParticleType::photon; }),
+    products_.erase(std::remove_if(products_.begin(), products_.end(),
+                      [](const auto& p) { return p.particle_.is_photon(); }),
       products_.end());
 
     // Determine product for D1S method
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
index 3ba2c0cfb05..ee560d6077e 100644
--- a/src/reaction_product.cpp
+++ b/src/reaction_product.cpp
@@ -26,7 +26,7 @@ ReactionProduct::ReactionProduct(hid_t group)
   // Read particle type
   std::string temp;
   read_attribute(group, "particle", temp);
-  particle_ = str_to_particle_type(temp);
+  particle_ = ParticleType {temp};
 
   // Read emission mode and decay rate
   read_attribute(group, "emission_mode", temp);
@@ -42,7 +42,7 @@ ReactionProduct::ReactionProduct(hid_t group)
   if (emission_mode_ == EmissionMode::delayed) {
     if (attribute_exists(group, "decay_rate")) {
       read_attribute(group, "decay_rate", decay_rate_);
-    } else if (particle_ == ParticleType::neutron) {
+    } else if (particle_.is_neutron()) {
       warning(fmt::format("Decay rate doesn't exist for delayed neutron "
                           "emission ({}).",
         object_name(group)));
@@ -85,7 +85,7 @@ ReactionProduct::ReactionProduct(hid_t group)
 
 ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)
 {
-  particle_ = ParticleType::photon;
+  particle_ = ParticleType::photon();
   emission_mode_ = EmissionMode::delayed;
 
   // Get chain nuclide object for radionuclide
diff --git a/src/simulation.cpp b/src/simulation.cpp
index b536ae5881f..18e40a8bc73 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -687,7 +687,7 @@ void initialize_data()
 
   for (const auto& nuc : data::nuclides) {
     if (nuc->grid_.size() >= 1) {
-      int neutron = static_cast<int>(ParticleType::neutron);
+      int neutron = ParticleType::neutron().transport_index();
       data::energy_min[neutron] =
         std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
       data::energy_max[neutron] =
@@ -698,7 +698,7 @@ void initialize_data()
   if (settings::photon_transport) {
     for (const auto& elem : data::elements) {
       if (elem->energy_.size() >= 1) {
-        int photon = static_cast<int>(ParticleType::photon);
+        int photon = ParticleType::photon().transport_index();
         int n = elem->energy_.size();
         data::energy_min[photon] =
           std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
@@ -711,9 +711,9 @@ void initialize_data()
       // Determine if minimum/maximum energy for bremsstrahlung is greater/less
       // than the current minimum/maximum
       if (data::ttb_e_grid.size() >= 1) {
-        int photon = static_cast<int>(ParticleType::photon);
-        int electron = static_cast<int>(ParticleType::electron);
-        int positron = static_cast<int>(ParticleType::positron);
+        int photon = ParticleType::photon().transport_index();
+        int electron = ParticleType::electron().transport_index();
+        int positron = ParticleType::positron().transport_index();
         int n_e = data::ttb_e_grid.size();
 
         const std::vector<int> charged = {electron, positron};
@@ -737,7 +737,7 @@ void initialize_data()
     // grid has not been allocated
     if (nuc->grid_.size() > 0) {
       double max_E = nuc->grid_[0].energy.back();
-      int neutron = static_cast<int>(ParticleType::neutron);
+      int neutron = ParticleType::neutron().transport_index();
       if (max_E == data::energy_max[neutron]) {
         write_message(7, "Maximum neutron transport energy: {} eV for {}",
           data::energy_max[neutron], nuc->name_);
@@ -754,7 +754,7 @@ void initialize_data()
   for (auto& nuc : data::nuclides) {
     nuc->init_grid();
   }
-  int neutron = static_cast<int>(ParticleType::neutron);
+  int neutron = ParticleType::neutron().transport_index();
   simulation::log_spacing =
     std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
     settings::n_log_bins;
diff --git a/src/source.cpp b/src/source.cpp
index b30b964d39c..8bd9c789352 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -37,6 +37,20 @@
 
 namespace openmc {
 
+namespace {
+
+void validate_particle_type(ParticleType type, const std::string& context)
+{
+  if (type.is_transportable())
+    return;
+
+  fatal_error(
+    fmt::format("Unsupported source particle type '{}' (PDG {}) in {}.",
+      type.str(), type.pdg_number(), context));
+}
+
+} // namespace
+
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -284,22 +298,15 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
 {
   // Check for particle type
   if (check_for_node(node, "particle")) {
-    auto temp_str = get_node_value(node, "particle", true, true);
-    if (temp_str == "neutron") {
-      particle_ = ParticleType::neutron;
-    } else if (temp_str == "photon") {
-      particle_ = ParticleType::photon;
+    auto temp_str = get_node_value(node, "particle", false, true);
+    particle_ = ParticleType(temp_str);
+    if (particle_ == ParticleType::photon() ||
+        particle_ == ParticleType::electron() ||
+        particle_ == ParticleType::positron()) {
       settings::photon_transport = true;
-    } else if (temp_str == "electron") {
-      particle_ = ParticleType::electron;
-      settings::photon_transport = true;
-    } else if (temp_str == "positron") {
-      particle_ = ParticleType::positron;
-      settings::photon_transport = true;
-    } else {
-      fatal_error(std::string("Unknown source particle type: ") + temp_str);
     }
   }
+  validate_particle_type(particle_, "IndependentSource");
 
   // Check for external source file
   if (check_for_node(node, "file")) {
@@ -390,7 +397,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
   // Sample energy and time for neutron and photon sources
   if (settings::solver_type != SolverType::RANDOM_RAY) {
     // Check for monoenergetic source above maximum particle energy
-    auto p = static_cast<int>(particle_);
+    auto p = particle_.transport_index();
     auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
     if (energy_ptr) {
       auto energies = xt::adapt(energy_ptr->x());
@@ -472,6 +479,11 @@ void FileSource::load_sites_from_file(const std::string& path)
     // Close file
     file_close(file_id);
   }
+
+  // Make sure particles in source file have valid types
+  for (const auto& site : this->sites_) {
+    validate_particle_type(site.particle, "FileSource");
+  }
 }
 
 SourceSite FileSource::sample(uint64_t* seed) const
@@ -585,6 +597,11 @@ MeshSource::MeshSource(pugi::xml_node node) : Source(node)
       std::make_unique<MeshElementSpatial>(mesh_idx, elem_index));
   }
 
+  // Make sure sources use valid particle types
+  for (const auto& src : sources_) {
+    validate_particle_type(src->particle_type(), "MeshSource");
+  }
+
   // the number of source distributions should either be one or equal to the
   // number of mesh elements
   if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) {
diff --git a/src/state_point.cpp b/src/state_point.cpp
index 8ccebeb05a7..455299529b3 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -22,6 +22,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
+#include "openmc/particle_type.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/tallies/derivative.h"
@@ -632,6 +633,7 @@ void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
   if (mpi::master || parallel) {
     file_id = file_open(filename_.c_str(), 'w', true);
     write_attribute(file_id, "filetype", "source");
+    write_attribute(file_id, "version", VERSION_STATEPOINT);
   }
 
   // Get pointer to source bank and write to file
@@ -677,6 +679,16 @@ std::string dtype_member_names(hid_t dtype_id)
 void read_source_bank(
   hid_t group_id, vector<SourceSite>& sites, bool distribute)
 {
+  bool legacy_particle_codes = true;
+  if (attribute_exists(group_id, "version")) {
+    array<int, 2> version;
+    read_attribute(group_id, "version", version);
+    if (version[0] > VERSION_STATEPOINT[0] ||
+        (version[0] == VERSION_STATEPOINT[0] && version[1] >= 2)) {
+      legacy_particle_codes = false;
+    }
+  }
+
   hid_t banktype = h5banktype(true);
 
   // Open the dataset
@@ -738,6 +750,12 @@ void read_source_bank(
     H5Sclose(memspace);
   H5Dclose(dset);
   H5Tclose(banktype);
+
+  if (legacy_particle_codes) {
+    for (auto& site : sites) {
+      site.particle = legacy_particle_index_to_type(site.particle.pdg_number());
+    }
+  }
 }
 
 void write_unstructured_mesh_results()
diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp
index eef1d1e63c5..031068f3a8d 100644
--- a/src/tallies/filter_particle.cpp
+++ b/src/tallies/filter_particle.cpp
@@ -13,7 +13,7 @@ void ParticleFilter::from_xml(pugi::xml_node node)
   // Convert to vector of ParticleType
   vector<ParticleType> types;
   for (auto& p : particles) {
-    types.push_back(str_to_particle_type(p));
+    types.emplace_back(p);
   }
   this->set_particles(types);
 }
@@ -47,7 +47,7 @@ void ParticleFilter::to_statepoint(hid_t filter_group) const
   Filter::to_statepoint(filter_group);
   vector<std::string> particles;
   for (auto p : particles_) {
-    particles.push_back(particle_type_to_str(p));
+    particles.push_back(p.str());
   }
   write_dataset(filter_group, "bins", particles);
 }
@@ -55,10 +55,10 @@ void ParticleFilter::to_statepoint(hid_t filter_group) const
 std::string ParticleFilter::text_label(int bin) const
 {
   const auto& p = particles_.at(bin);
-  return fmt::format("Particle: {}", particle_type_to_str(p));
+  return fmt::format("Particle: {}", p.str());
 }
 
-extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[])
+extern "C" int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[])
 {
   if (int err = verify_filter(idx))
     return err;
@@ -68,7 +68,7 @@ extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[])
   if (pf) {
     const auto& particles = pf->particles();
     for (int i = 0; i < particles.size(); i++) {
-      bins[i] = static_cast<int>(particles[i]);
+      bins[i] = particles[i].pdg_number();
     }
   } else {
     set_errmsg("The filter at the specified index is not a ParticleFilter");
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 6eef1da9cfd..2432f5c2f18 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -289,12 +289,12 @@ Tally::Tally(pugi::xml_node node)
       const auto& f = model::tally_filters[particle_filter_index].get();
       auto pf = dynamic_cast<ParticleFilter*>(f);
       for (auto p : pf->particles()) {
-        if (p != ParticleType::neutron) {
+        if (!p.is_neutron()) {
           warning(fmt::format(
             "Particle filter other than NEUTRON used with "
             "photon transport turned off. All tallies for particle type {}"
             " will have no scores",
-            static_cast<int>(p)));
+            p.str()));
         }
       }
     }
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 67e851644a9..51b5d9ffcc5 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -328,10 +328,10 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux,
 //! Helper function to obtain reaction Q value for photons and charged particles
 double get_reaction_q_value(const Particle& p)
 {
-  if (p.type() == ParticleType::photon && p.event_mt() == PAIR_PROD) {
+  if (p.type().is_photon() && p.event_mt() == PAIR_PROD) {
     // pair production
     return -2 * MASS_ELECTRON_EV;
-  } else if (p.type() == ParticleType::positron) {
+  } else if (p.type() == ParticleType::positron()) {
     // positron annihilation
     return 2 * MASS_ELECTRON_EV;
   } else {
@@ -344,7 +344,7 @@ double get_reaction_q_value(const Particle& p)
 double score_particle_heating(const Particle& p, const Tally& tally,
   double flux, int rxn_bin, int i_nuclide, double atom_density)
 {
-  if (p.type() == ParticleType::neutron)
+  if (p.type().is_neutron())
     return score_neutron_heating(
       p, tally, flux, rxn_bin, i_nuclide, atom_density);
   if (i_nuclide == -1 || i_nuclide == p.event_nuclide() ||
@@ -584,8 +584,6 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
   // Get the pre-collision energy of the particle.
   auto E = p.E_last();
 
-  using Type = ParticleType;
-
   for (auto i = 0; i < tally.scores_.size(); ++i) {
     auto score_bin = tally.scores_[i];
     auto score_index = start_index + i;
@@ -598,9 +596,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 
     case SCORE_TOTAL:
       if (i_nuclide >= 0) {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = p.neutron_xs(i_nuclide).total * atom_density * flux;
-        } else if (p.type() == Type::photon) {
+        } else if (p.type().is_photon()) {
           score = p.photon_xs(i_nuclide).total * atom_density * flux;
         }
       } else {
@@ -609,7 +607,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_INVERSE_VELOCITY:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Score inverse velocity in units of s/cm.
@@ -617,11 +615,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_SCATTER:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       if (i_nuclide >= 0) {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           const auto& micro = p.neutron_xs(i_nuclide);
           score = (micro.total - micro.absorption) * atom_density * flux;
         } else {
@@ -629,7 +627,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
           score = (micro.coherent + micro.incoherent) * atom_density * flux;
         }
       } else {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = (p.macro_xs().total - p.macro_xs().absorption) * flux;
         } else {
           score = (p.macro_xs().coherent + p.macro_xs().incoherent) * flux;
@@ -638,11 +636,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_ABSORPTION:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       if (i_nuclide >= 0) {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = p.neutron_xs(i_nuclide).absorption * atom_density * flux;
         } else {
           const auto& xs = p.photon_xs(i_nuclide);
@@ -650,7 +648,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
             (xs.total - xs.coherent - xs.incoherent) * atom_density * flux;
         }
       } else {
-        if (p.type() == Type::neutron) {
+        if (p.type().is_neutron()) {
           score = p.macro_xs().absorption * flux;
         } else {
           score =
@@ -806,7 +804,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
           // this loop.
           for (auto d = 1; d < rxn.products_.size(); ++d) {
             const auto& product = rxn.products_[d];
-            if (product.particle_ != Type::neutron)
+            if (!product.particle_.is_neutron())
               continue;
 
             auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
@@ -860,7 +858,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
                 // this loop.
                 for (auto d = 1; d < rxn.products_.size(); ++d) {
                   const auto& product = rxn.products_[d];
-                  if (product.particle_ != Type::neutron)
+                  if (!product.particle_.is_neutron())
                     continue;
 
                   auto yield =
@@ -911,7 +909,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       continue;
 
     case ELASTIC:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       if (i_nuclide >= 0) {
@@ -943,7 +941,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 
     case SCORE_IFP_TIME_NUM:
       if (settings::ifp_on) {
-        if ((p.type() == Type::neutron) && (p.fission())) {
+        if (p.type().is_neutron() && p.fission()) {
           if (is_generation_time_or_both()) {
             const auto& lifetimes =
               simulation::ifp_source_lifetime_bank[p.current_work() - 1];
@@ -957,7 +955,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 
     case SCORE_IFP_BETA_NUM:
       if (settings::ifp_on) {
-        if ((p.type() == Type::neutron) && (p.fission())) {
+        if (p.type().is_neutron() && p.fission()) {
           if (is_beta_effective_or_both()) {
             const auto& delayed_groups =
               simulation::ifp_source_delayed_group_bank[p.current_work() - 1];
@@ -982,7 +980,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 
     case SCORE_IFP_DENOM:
       if (settings::ifp_on) {
-        if ((p.type() == Type::neutron) && (p.fission())) {
+        if (p.type().is_neutron() && p.fission()) {
           int ifp_data_size;
           if (is_beta_effective_or_both()) {
             ifp_data_size = static_cast<int>(
@@ -1012,7 +1010,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       if (!simulation::need_depletion_rx)
         goto default_case;
 
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       int m;
@@ -1045,7 +1043,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
     case INCOHERENT:
     case PHOTOELECTRIC:
     case PAIR_PROD:
-      if (p.type() != Type::photon)
+      if (!p.type().is_photon())
         continue;
 
       if (i_nuclide >= 0) {
@@ -1075,7 +1073,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 
       // The default block is really only meant for redundant neutron reactions
       // (e.g. 444, 901)
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Any other cross section has to be calculated on-the-fly
@@ -1129,8 +1127,6 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
                             p.neutron_xs(p.event_nuclide()).total
                         : 0.0;
 
-  using Type = ParticleType;
-
   for (auto i = 0; i < tally.scores_.size(); ++i) {
     auto score_bin = tally.scores_[i];
     auto score_index = start_index + i;
@@ -1141,7 +1137,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       // All events score to a flux bin. We actually use a collision estimator
       // in place of an analog one since there is no way to count 'events'
       // exactly for the flux
-      if (p.type() == Type::neutron || p.type() == Type::photon) {
+      if (p.type().is_neutron() || p.type().is_photon()) {
         score = flux * p.wgt_last() / p.macro_xs().total;
       } else {
         score = 0.;
@@ -1155,7 +1151,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_INVERSE_VELOCITY:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // All events score to an inverse velocity bin. We actually use a
@@ -1165,7 +1161,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_SCATTER:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       // Skip any event where the particle didn't scatter
@@ -1177,7 +1173,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_NU_SCATTER:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Only analog estimators are available.
@@ -1202,7 +1198,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case SCORE_ABSORPTION:
-      if (p.type() != Type::neutron && p.type() != Type::photon)
+      if (!p.type().is_neutron() && !p.type().is_photon())
         continue;
 
       if (settings::survival_biasing) {
@@ -1431,7 +1427,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
             // this loop.
             for (auto d = 1; d < rxn.products_.size(); ++d) {
               const auto& product = rxn.products_[d];
-              if (product.particle_ != Type::neutron)
+              if (!product.particle_.is_neutron())
                 continue;
 
               auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
@@ -1523,7 +1519,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       continue;
 
     case ELASTIC:
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Check if event MT matches
@@ -1552,7 +1548,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       if (!simulation::need_depletion_rx)
         goto default_case;
 
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Check if the event MT matches
@@ -1565,7 +1561,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
     case INCOHERENT:
     case PHOTOELECTRIC:
     case PAIR_PROD:
-      if (p.type() != Type::photon)
+      if (!p.type().is_photon())
         continue;
 
       if (score_bin == PHOTOELECTRIC) {
@@ -1592,7 +1588,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
 
       // The default block is really only meant for redundant neutron reactions
       // (e.g. 444, 901)
-      if (p.type() != Type::neutron)
+      if (!p.type().is_neutron())
         continue;
 
       // Any other score is assumed to be a MT number. Thus, we just need
@@ -2316,10 +2312,7 @@ void score_analog_tally_ce(Particle& p)
   // Since electrons/positrons are not transported, we assign a flux of zero.
   // Note that the heating score does NOT use the flux and will be non-zero for
   // electrons/positrons.
-  double flux =
-    (p.type() == ParticleType::neutron || p.type() == ParticleType::photon)
-      ? 1.0
-      : 0.0;
+  double flux = (p.type().is_neutron() || p.type().is_photon()) ? 1.0 : 0.0;
 
   for (auto i_tally : model::active_analog_tallies) {
     const Tally& tally {*model::tallies[i_tally]};
@@ -2447,7 +2440,7 @@ void score_tracklength_tally_general(
             if (j == C_NONE) {
               // Determine log union grid index
               if (i_log_union == C_NONE) {
-                int neutron = static_cast<int>(ParticleType::neutron);
+                int neutron = ParticleType::neutron().transport_index();
                 i_log_union = std::log(p.E() / data::energy_min[neutron]) /
                               simulation::log_spacing;
               }
@@ -2544,7 +2537,7 @@ void score_collision_tally(Particle& p)
 {
   // Determine the collision estimate of the flux
   double flux = 0.0;
-  if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) {
+  if (p.type().is_neutron() || p.type().is_photon()) {
     flux = p.wgt_last() / p.macro_xs().total;
   }
 
@@ -2578,7 +2571,7 @@ void score_collision_tally(Particle& p)
           if (j == C_NONE) {
             // Determine log union grid index
             if (i_log_union == C_NONE) {
-              int neutron = static_cast<int>(ParticleType::neutron);
+              int neutron = ParticleType::neutron().transport_index();
               i_log_union = std::log(p.E() / data::energy_min[neutron]) /
                             simulation::log_spacing;
             }
diff --git a/src/track_output.cpp b/src/track_output.cpp
index f4344d50f74..e86f774f046 100644
--- a/src/track_output.cpp
+++ b/src/track_output.cpp
@@ -122,7 +122,7 @@ void finalize_particle_track(Particle& p)
   int offset = 0;
   for (auto& track_i : p.tracks()) {
     offsets.push_back(offset);
-    particles.push_back(static_cast<int>(track_i.particle));
+    particles.push_back(track_i.particle.pdg_number());
     offset += track_i.states.size();
     tracks.insert(tracks.end(), track_i.states.begin(), track_i.states.end());
   }
diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp
index 4838e459152..5ca4addbbf7 100644
--- a/src/weight_windows.cpp
+++ b/src/weight_windows.cpp
@@ -59,7 +59,7 @@ void apply_weight_windows(Particle& p)
     return;
 
   // WW on photon and neutron only
-  if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon)
+  if (!p.type().is_neutron() && !p.type().is_photon())
     return;
 
   // skip dead or no energy
@@ -179,7 +179,7 @@ WeightWindows::WeightWindows(pugi::xml_node node)
 
   // get the particle type
   auto particle_type_str = std::string(get_node_value(node, "particle_type"));
-  particle_type_ = openmc::str_to_particle_type(particle_type_str);
+  particle_type_ = ParticleType {particle_type_str};
 
   // Determine associated mesh
   int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
@@ -252,7 +252,7 @@ WeightWindows* WeightWindows::from_hdf5(
 
   std::string particle_type;
   read_dataset(ww_group, "particle_type", particle_type);
-  wws->particle_type_ = openmc::str_to_particle_type(particle_type);
+  wws->particle_type_ = ParticleType {particle_type};
 
   read_dataset<double>(ww_group, "energy_bounds", wws->energy_bounds_);
 
@@ -284,7 +284,10 @@ void WeightWindows::set_defaults()
 {
   // set energy bounds to the min/max energy supported by the data
   if (energy_bounds_.size() == 0) {
-    int p_type = static_cast<int>(particle_type_);
+    int p_type = particle_type_.transport_index();
+    if (p_type == C_NONE) {
+      fatal_error("Weight windows particle is not supported for transport.");
+    }
     energy_bounds_.push_back(data::energy_min[p_type]);
     energy_bounds_.push_back(data::energy_max[p_type]);
   }
@@ -345,10 +348,9 @@ void WeightWindows::set_energy_bounds(span<const double> bounds)
 
 void WeightWindows::set_particle_type(ParticleType p_type)
 {
-  if (p_type != ParticleType::neutron && p_type != ParticleType::photon)
-    fatal_error(
-      fmt::format("Particle type '{}' cannot be applied to weight windows.",
-        particle_type_to_str(p_type)));
+  if (!p_type.is_neutron() && !p_type.is_photon())
+    fatal_error(fmt::format(
+      "Particle type '{}' cannot be applied to weight windows.", p_type.str()));
   particle_type_ = p_type;
 }
 
@@ -608,8 +610,7 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value,
     if (p_it == particles.end()) {
       auto msg = fmt::format("Particle type '{}' not present on Filter {} for "
                              "Tally {} used to update WeightWindows {}",
-        particle_type_to_str(this->particle_type_), pf->id(), tally->id(),
-        this->id());
+        this->particle_type_.str(), pf->id(), tally->id(), this->id());
       fatal_error(msg);
     }
 
@@ -818,8 +819,7 @@ void WeightWindows::to_hdf5(hid_t group) const
   hid_t ww_group = create_group(group, fmt::format("weight_windows_{}", id()));
 
   write_dataset(ww_group, "mesh", this->mesh()->id());
-  write_dataset(
-    ww_group, "particle_type", openmc::particle_type_to_str(particle_type_));
+  write_dataset(ww_group, "particle_type", particle_type_.str());
   write_dataset(ww_group, "energy_bounds", energy_bounds_);
   write_dataset(ww_group, "lower_ww_bounds", lower_ww_);
   write_dataset(ww_group, "upper_ww_bounds", upper_ww_);
@@ -846,8 +846,8 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
         max_realizations_, active_batches);
     warning(msg);
   }
-  auto tmp_str = get_node_value(node, "particle_type", true, true);
-  auto particle_type = str_to_particle_type(tmp_str);
+  auto tmp_str = get_node_value(node, "particle_type", false, true);
+  auto particle_type = ParticleType {tmp_str};
 
   update_interval_ = std::stoi(get_node_value(node, "update_interval"));
   on_the_fly_ = get_node_value_bool(node, "on_the_fly");
@@ -856,7 +856,10 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
   if (check_for_node(node, "energy_bounds")) {
     e_bounds = get_node_array<double>(node, "energy_bounds");
   } else {
-    int p_type = static_cast<int>(particle_type);
+    int p_type = particle_type.transport_index();
+    if (p_type == C_NONE) {
+      fatal_error("Weight windows particle is not supported for transport.");
+    }
     e_bounds.push_back(data::energy_min[p_type]);
     e_bounds.push_back(data::energy_max[p_type]);
   }
@@ -1110,23 +1113,25 @@ extern "C" int openmc_weight_windows_get_energy_bounds(
   return 0;
 }
 
-extern "C" int openmc_weight_windows_set_particle(int32_t index, int particle)
+extern "C" int openmc_weight_windows_set_particle(
+  int32_t index, int32_t particle)
 {
   if (int err = verify_ww_index(index))
     return err;
 
   const auto& wws = variance_reduction::weight_windows.at(index);
-  wws->set_particle_type(static_cast<ParticleType>(particle));
+  wws->set_particle_type(ParticleType {particle});
   return 0;
 }
 
-extern "C" int openmc_weight_windows_get_particle(int32_t index, int* particle)
+extern "C" int openmc_weight_windows_get_particle(
+  int32_t index, int32_t* particle)
 {
   if (int err = verify_ww_index(index))
     return err;
 
   const auto& wws = variance_reduction::weight_windows.at(index);
-  *particle = static_cast<int>(wws->particle_type());
+  *particle = wws->particle_type().pdg_number();
   return 0;
 }
 
diff --git a/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp b/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp
index 909830e035c..d0c26f26c2a 100644
--- a/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp
+++ b/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp
@@ -25,7 +25,7 @@ TEST_CASE("MCPL stat:sum field")
 
     // Initialize test particles
     for (int i = 0; i < 100; ++i) {
-      source_bank[i].particle = openmc::ParticleType::neutron;
+      source_bank[i].particle = openmc::ParticleType::neutron();
       source_bank[i].r = {i * 0.1, i * 0.2, i * 0.3};
       source_bank[i].u = {0.0, 0.0, 1.0};
       source_bank[i].E = 2.0e6; // 2 MeV
@@ -63,7 +63,7 @@ TEST_CASE("MCPL stat:sum field")
 
       // Initialize particles
       for (int i = 0; i < count; ++i) {
-        source_bank[i].particle = openmc::ParticleType::neutron;
+        source_bank[i].particle = openmc::ParticleType::neutron();
         source_bank[i].r = {0.0, 0.0, 0.0};
         source_bank[i].u = {0.0, 0.0, 1.0};
         source_bank[i].E = 1.0e6;
diff --git a/tests/regression_tests/source_dlopen/source_sampling.cpp b/tests/regression_tests/source_dlopen/source_sampling.cpp
index fc61ef1fdd8..678eea4be06 100644
--- a/tests/regression_tests/source_dlopen/source_sampling.cpp
+++ b/tests/regression_tests/source_dlopen/source_sampling.cpp
@@ -10,7 +10,7 @@ class CustomSource : public openmc::Source {
   {
     openmc::SourceSite particle;
     // wgt
-    particle.particle = openmc::ParticleType::neutron;
+    particle.particle = openmc::ParticleType::neutron();
     particle.wgt = 1.0;
     // position
 
diff --git a/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp b/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp
index bf49af4c388..81f3770c62b 100644
--- a/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp
+++ b/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp
@@ -2,36 +2,37 @@
 #include "openmc/source.h"
 
 class CustomSource : public openmc::Source {
-  public:
-    CustomSource(double energy) : energy_(energy) { }
+public:
+  CustomSource(double energy) : energy_(energy) {}
 
-    // Samples from an instance of this class.
-    openmc::SourceSite sample(uint64_t* seed) const
-    {
-      openmc::SourceSite particle;
-      // wgt
-      particle.particle = openmc::ParticleType::neutron;
-      particle.wgt = 1.0;
-      // position
-      particle.r.x = 0.0;
-      particle.r.y = 0.0;
-      particle.r.z = 0.0;
-      // angle
-      particle.u = {1.0, 0.0, 0.0};
-      particle.E = this->energy_;
-      particle.delayed_group = 0;
+  // Samples from an instance of this class.
+  openmc::SourceSite sample(uint64_t* seed) const
+  {
+    openmc::SourceSite particle;
+    // wgt
+    particle.particle = openmc::ParticleType::neutron();
+    particle.wgt = 1.0;
+    // position
+    particle.r.x = 0.0;
+    particle.r.y = 0.0;
+    particle.r.z = 0.0;
+    // angle
+    particle.u = {1.0, 0.0, 0.0};
+    particle.E = this->energy_;
+    particle.delayed_group = 0;
 
-      return particle;
-    }
+    return particle;
+  }
 
-  private:
-    double energy_;
+private:
+  double energy_;
 };
 
-// A function to create a unique pointer to an instance of this class when generated
-// via a plugin call using dlopen/dlsym.
-// You must have external C linkage here otherwise dlopen will not find the file
-extern "C" std::unique_ptr<CustomSource> openmc_create_source(std::string parameter)
+// A function to create a unique pointer to an instance of this class when
+// generated via a plugin call using dlopen/dlsym. You must have external C
+// linkage here otherwise dlopen will not find the file
+extern "C" std::unique_ptr<CustomSource> openmc_create_source(
+  std::string parameter)
 {
   double energy = std::stod(parameter);
   return std::make_unique<CustomSource>(energy);
diff --git a/tests/regression_tests/surface_source/surface_source_true.h5 b/tests/regression_tests/surface_source/surface_source_true.h5
index d343d12aed6..02e94c90b99 100644
Binary files a/tests/regression_tests/surface_source/surface_source_true.h5 and b/tests/regression_tests/surface_source/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/test.py b/tests/regression_tests/surface_source_write/test.py
index f144eb82a73..094df1f8b84 100644
--- a/tests/regression_tests/surface_source_write/test.py
+++ b/tests/regression_tests/surface_source_write/test.py
@@ -634,7 +634,7 @@ def return_surface_source_data(filepath):
         wgt = point.wgt
         delayed_group = point.delayed_group
         surf_id = point.surf_id
-        particle = point.particle
+        particle = point.particle.pdg_number
         key = (
             f"{r[0]:.10e} {r[1]:.10e} {r[2]:.10e} {u[0]:.10e} {u[1]:.10e} {u[2]:.10e}"
             f"{e:.10e} {time:.10e} {wgt:.10e} {delayed_group} {surf_id} {particle}"
diff --git a/tests/unit_tests/test_particle_type.py b/tests/unit_tests/test_particle_type.py
new file mode 100644
index 00000000000..67e385ea632
--- /dev/null
+++ b/tests/unit_tests/test_particle_type.py
@@ -0,0 +1,298 @@
+"""Unit tests for ParticleType class."""
+
+import pytest
+from openmc import ParticleType
+
+
+# Tests for creating ParticleType instances
+
+def test_create_from_int():
+    """Test creation from PDG number."""
+    p = ParticleType(2112)
+    assert p.pdg_number == 2112
+    assert int(p) == 2112
+
+
+def test_create_from_string_name():
+    """Test creation from particle name."""
+    p = ParticleType('neutron')
+    assert p.pdg_number == 2112
+
+    p = ParticleType('photon')
+    assert p.pdg_number == 22
+
+
+def test_create_from_string_aliases():
+    """Test creation from particle aliases."""
+    assert ParticleType('n').pdg_number == 2112
+    assert ParticleType('gamma').pdg_number == 22
+    assert ParticleType('p').pdg_number == 2212
+    assert ParticleType('proton').pdg_number == 2212
+    assert ParticleType('d').pdg_number == 1000010020
+    assert ParticleType('t').pdg_number == 1000010030
+
+
+def test_create_from_string_nuclide():
+    """Test creation from GNDS nuclide name."""
+    p = ParticleType('He4')
+    assert p.pdg_number == 1000020040
+
+    p = ParticleType('U235')
+    assert p.pdg_number == 1000922350
+
+    p = ParticleType('Am242_m1')
+    assert p.pdg_number == 1000952421
+
+
+def test_create_from_string_pdg_prefix():
+    """Test creation with pdg: prefix."""
+    p = ParticleType('pdg:2112')
+    assert p.pdg_number == 2112
+
+    p = ParticleType('PDG:22')
+    assert p.pdg_number == 22
+
+
+def test_create_from_particle_type():
+    """Test creation from existing ParticleType."""
+    p1 = ParticleType(2112)
+    p2 = ParticleType(p1)
+    assert p1 == p2
+    assert p1 is not p2  # Different instances
+
+
+def test_legacy_particle_indices():
+    """Test backward compatibility with legacy indices 0-3."""
+    assert ParticleType(0) == ParticleType.NEUTRON
+    assert ParticleType(1) == ParticleType.PHOTON
+    assert ParticleType(2) == ParticleType.ELECTRON
+    assert ParticleType(3) == ParticleType.POSITRON
+
+
+def test_create_invalid_type():
+    """Test creation with invalid type raises TypeError."""
+    with pytest.raises(TypeError):
+        ParticleType([2112])
+
+    with pytest.raises(TypeError):
+        ParticleType({'pdg': 2112})
+
+
+def test_create_invalid_string():
+    """Test creation with invalid string raises ValueError."""
+    with pytest.raises(ValueError):
+        ParticleType('')
+
+    with pytest.raises(ValueError):
+        ParticleType('pdg:invalid')
+
+
+def test_create_case_insensitive():
+    """Test that string parsing is case insensitive for aliases."""
+    assert ParticleType('NEUTRON').pdg_number == 2112
+    assert ParticleType('Neutron').pdg_number == 2112
+    assert ParticleType('PHOTON').pdg_number == 22
+
+
+# Tests for equality and comparison
+
+def test_equality_same_pdg():
+    """Test that instances with same PDG are equal."""
+    p1 = ParticleType(2112)
+    p2 = ParticleType(2112)
+    assert p1 == p2
+
+
+def test_equality_int():
+    """Test equality comparison with int."""
+    p = ParticleType(2112)
+    assert p == 2112
+    assert 2112 == p
+
+
+def test_equality_string():
+    """Test equality comparison with string."""
+    p = ParticleType(2112)
+    assert p == 'neutron'
+    assert p == 'pdg:2112'
+    assert p == 'n'
+
+
+def test_inequality():
+    """Test inequality comparisons."""
+    p1 = ParticleType(2112)
+    p2 = ParticleType(22)
+    assert p1 != p2
+    assert p1 != 22
+    assert p1 != 'photon'
+
+
+def test_equality_with_constants():
+    """Test equality with class constants."""
+    p = ParticleType(2112)
+    assert p == ParticleType.NEUTRON
+    assert ParticleType.NEUTRON == p
+
+
+def test_equality_invalid_string():
+    """Test equality with invalid string returns False."""
+    p = ParticleType(2112)
+    assert not (p == 'invalid_particle')
+    assert p != 'invalid_particle'
+
+
+# Tests for hashing behavior
+
+def test_hash_consistency():
+    """Test that equal instances hash to the same value."""
+    p1 = ParticleType(2112)
+    p2 = ParticleType(2112)
+    assert hash(p1) == hash(p2)
+
+
+def test_set_deduplication():
+    """Test that equal instances deduplicate in sets."""
+    p1 = ParticleType(2112)
+    p2 = ParticleType(2112)
+    s = {p1, p2}
+    assert len(s) == 1
+
+
+def test_dict_key():
+    """Test use as dictionary key."""
+    p1 = ParticleType(2112)
+    p2 = ParticleType(2112)
+    d = {p1: 'neutron'}
+    assert d[p2] == 'neutron'
+
+
+def test_hash_different_particles():
+    """Test that different particles have different hashes (usually)."""
+    p1 = ParticleType(2112)
+    p2 = ParticleType(22)
+    # Different PDG numbers should (almost always) have different hashes
+    assert hash(p1) != hash(p2)
+
+
+# Tests for properties and computed attributes
+
+def test_pdg_number_property():
+    """Test pdg_number property."""
+    p = ParticleType(2112)
+    assert p.pdg_number == 2112
+    assert isinstance(p.pdg_number, int)
+
+
+def test_int_conversion():
+    """Test __int__ conversion."""
+    p = ParticleType(1000020040)
+    assert int(p) == 1000020040
+
+
+def test_zam_elementary():
+    """Test zam property for elementary particles."""
+    assert ParticleType.NEUTRON.zam is None
+    assert ParticleType.PHOTON.zam is None
+    assert ParticleType.ELECTRON.zam is None
+    assert ParticleType.POSITRON.zam is None
+
+
+def test_zam_nucleus():
+    """Test zam property for nuclear particles."""
+    he4 = ParticleType('He4')
+    assert he4.zam == (2, 4, 0)
+
+    u235 = ParticleType('U235')
+    Z, A, m = u235.zam
+    assert Z == 92
+    assert A == 235
+    assert m == 0
+
+
+def test_zam_metastable():
+    """Test zam property for metastable nuclei."""
+    # Am242m has m=1
+    am242m = ParticleType('Am242_m1')
+    Z, A, m = am242m.zam
+    assert Z == 95
+    assert A == 242
+    assert m == 1
+
+
+def test_is_nucleus_false():
+    """Test is_nucleus for elementary particles."""
+    assert not ParticleType.NEUTRON.is_nucleus
+    assert not ParticleType.PHOTON.is_nucleus
+    assert not ParticleType.ELECTRON.is_nucleus
+    assert not ParticleType.POSITRON.is_nucleus
+
+
+def test_is_nucleus_true():
+    """Test is_nucleus for nuclear particles."""
+    assert ParticleType.ALPHA.is_nucleus
+    assert ParticleType('He4').is_nucleus
+    assert ParticleType('U235').is_nucleus
+    assert ParticleType.DEUTERON.is_nucleus
+    assert ParticleType.TRITON.is_nucleus
+
+
+# Tests for __str__ and __repr__
+
+def test_str_elementary():
+    """Test string representation of elementary particles."""
+    assert str(ParticleType.NEUTRON) == 'neutron'
+    assert str(ParticleType.PHOTON) == 'photon'
+    assert str(ParticleType.ELECTRON) == 'electron'
+    assert str(ParticleType.POSITRON) == 'positron'
+    assert str(ParticleType.PROTON) == 'H1'
+
+
+def test_str_nucleus():
+    """Test string representation of nuclei."""
+    assert str(ParticleType.ALPHA) == 'He4'
+    assert str(ParticleType('U235')) == 'U235'
+    assert str(ParticleType.DEUTERON) == 'H2'
+    assert str(ParticleType.TRITON) == 'H3'
+
+
+def test_str_arbitrary_pdg():
+    """Test string representation of arbitrary PDG number."""
+    # PDG number that doesn't match any known particle
+    p = ParticleType(12345)
+    assert str(p) == 'pdg:12345'
+
+
+def test_repr():
+    """Test repr includes PDG number."""
+    p = ParticleType.NEUTRON
+    repr_str = repr(p)
+    assert 'ParticleType' in repr_str
+    assert 'PDG=2112' in repr_str
+    assert 'neutron' in repr_str
+
+
+def test_repr_nucleus():
+    """Test repr for nuclear particles."""
+    p = ParticleType.ALPHA
+    repr_str = repr(p)
+    assert 'ParticleType' in repr_str
+    assert 'PDG=1000020040' in repr_str
+    assert 'He4' in repr_str
+
+
+def test_create_from_numpy_int():
+    """Test creation from numpy integer types."""
+    import numpy as np
+    p = ParticleType(np.int32(2112))
+    assert p.pdg_number == 2112
+
+    p = ParticleType(np.int64(22))
+    assert p.pdg_number == 22
+
+
+def test_equality_numpy_int():
+    """Test equality comparison with numpy integer types."""
+    import numpy as np
+    p = ParticleType(2112)
+    assert p == np.int32(2112)
+    assert p == np.int64(2112)
diff --git a/tests/unit_tests/test_source_file.py b/tests/unit_tests/test_source_file.py
index 41906c80f83..3e84fbe5176 100644
--- a/tests/unit_tests/test_source_file.py
+++ b/tests/unit_tests/test_source_file.py
@@ -42,7 +42,7 @@ def test_source_file(run_in_tmpdir):
     assert np.all(arr['E'] == n - np.arange(n))
     assert np.all(arr['wgt'] == 1.0)
     assert np.all(arr['delayed_group'] == 0)
-    assert np.all(arr['particle'] == 0)
+    assert np.all(arr['particle'] == 2112)  # PDG number for neutron
 
     # Ensure sites read in are consistent
     sites = openmc.ParticleList.from_hdf5('test_source.h5')
@@ -64,7 +64,7 @@ def test_source_file(run_in_tmpdir):
     dgs = np.array([s.delayed_group for s in sites])
     assert np.all(dgs == 0)
     p_types = np.array([s.particle for s in sites])
-    assert np.all(p_types == 0)
+    assert np.all(p_types == 2112)  # PDG number for neutron
 
     # Ensure a ParticleList item is a SourceParticle
     site = sites[0]
diff --git a/tests/unit_tests/test_tracks.py b/tests/unit_tests/test_tracks.py
index fa866415984..74da83552f6 100644
--- a/tests/unit_tests/test_tracks.py
+++ b/tests/unit_tests/test_tracks.py
@@ -70,7 +70,7 @@ def test_tracks(sphere_model, particle, run_in_tmpdir):
         particle_track = track.particle_tracks[0]
         assert isinstance(particle_track, openmc.ParticleTrack)
 
-        assert particle_track.particle.name.lower() == particle
+        assert str(particle_track.particle).lower() == particle
         assert isinstance(particle_track.states, np.ndarray)
 
         # Sanity checks on actual data
@@ -140,8 +140,10 @@ def test_filter(sphere_model, run_in_tmpdir):
     assert matches == tracks
     matches = tracks.filter(state_filter=lambda s: s['E'] > 0.0)
     assert matches == tracks
-    matches = tracks.filter(particle='bunnytron')
+    matches = tracks.filter(particle='proton')
     assert matches == []
+    with pytest.raises(ValueError):
+        tracks.filter(particle='bunnytron')
 
 
 def test_write_to_vtk(sphere_model):
diff --git a/tests/unit_tests/weightwindows/test_wwinp_reader.py b/tests/unit_tests/weightwindows/test_wwinp_reader.py
index 637463bb2e6..601517d47ca 100644
--- a/tests/unit_tests/weightwindows/test_wwinp_reader.py
+++ b/tests/unit_tests/weightwindows/test_wwinp_reader.py
@@ -52,7 +52,7 @@
 n_e_bounds = (np.array([0.0,
                          100000.0,
                          146780.0]),)
-n_particles = ('neutron',)
+n_particles = [openmc.ParticleType.NEUTRON]
 
 # expected results - neutron and photon data
 np_mesh = openmc.RectilinearMesh()
@@ -63,7 +63,7 @@
 
 np_e_bounds = (np.array([0.0, 100000.0, 146780.0, 215440.0]),
                np.array([0.0, 1.0E8]))
-np_particles = ('neutron', 'photon')
+np_particles = [openmc.ParticleType.NEUTRON, openmc.ParticleType.PHOTON]
 
 # expected results - photon data only
 p_mesh = openmc.RectilinearMesh()
@@ -74,7 +74,7 @@
 p_mesh.z_grid = np.array([-50.0, 50.0])
 
 p_e_bounds = (np.array([0.0, 100000.0, 146780.0, 215440.0, 316230.0]),)
-p_particles = ('photon',)
+p_particles = [openmc.ParticleType.PHOTON]
 
 expected_results = [('wwinp_n', n_mesh, n_particles, n_e_bounds),
                     ('wwinp_np', np_mesh, np_particles, np_e_bounds),