diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 68abbd05ec..b0ebd0a915 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -17,6 +17,8 @@ Next release
 
 *Fixed*
 
+* Divide by zero with MPCD if rigid body has constituent with zero mass (#2228)
+
 6.1.0 (2026-02-07)
 ^^^^^^^^^^^^^^^^^^^^
 
diff --git a/hoomd/mpcd/CollisionMethod.cc b/hoomd/mpcd/CollisionMethod.cc
index 23c7d5c6e6..9e86191e56 100644
--- a/hoomd/mpcd/CollisionMethod.cc
+++ b/hoomd/mpcd/CollisionMethod.cc
@@ -437,8 +437,14 @@ void mpcd::CollisionMethod::beginThermalizeConstituentParticles(uint64_t timeste
             }
 
         const Scalar mass_const = h_velocity.data[idx].w;
-        const unsigned int tag = h_tag.data[idx];
+        // don't thermalize particles with zero mass
+        if (mass_const == Scalar(0))
+            {
+            continue;
+            }
+
         // draw random velocities from normal distribution
+        const unsigned int tag = h_tag.data[idx];
         hoomd::RandomGenerator rng(
             hoomd::Seed(hoomd::RNGIdentifier::CollisionMethod, timestep, seed),
             hoomd::Counter(tag, 1));
@@ -571,6 +577,11 @@ void mpcd::CollisionMethod::finishThermalizeConstituentParticles(uint64_t timest
 
         // get velocities and masses
         Scalar4 vel_constituent = h_velocity.data[idx];
+        if (vel_constituent.w == Scalar(0))
+            {
+            // constituents with zero mass don't get thermalized
+            continue;
+            }
         const Scalar4 thermal_vel_mass = h_alt_vel.data[idx];
         vec3<Scalar> thermal_vel(thermal_vel_mass);
 
diff --git a/hoomd/mpcd/CollisionMethod.cu b/hoomd/mpcd/CollisionMethod.cu
index b476476cb4..73681636da 100644
--- a/hoomd/mpcd/CollisionMethod.cu
+++ b/hoomd/mpcd/CollisionMethod.cu
@@ -47,6 +47,11 @@ __global__ void draw_velocities_constituent_particles(Scalar3* d_linmom_accum,
         }
 
     const Scalar mass_const = d_velocity[idx].w;
+    // don't thermalize particles with zero mass
+    if (mass_const == Scalar(0))
+        {
+        return;
+        }
     const unsigned int tag = d_tag[idx];
     // draw random velocities from normal distribution
     hoomd::RandomGenerator rng(hoomd::Seed(hoomd::RNGIdentifier::CollisionMethod, timestep, seed),
@@ -167,6 +172,11 @@ __global__ void apply_thermalized_velocity_vectors(const Scalar3* d_angmom_accum
 
     // get velocities and masses
     Scalar4 vel_constituent = d_velocity[idx];
+    if (vel_constituent.w == Scalar(0))
+        {
+        // constituents with zero mass don't get thermalized
+        return;
+        }
     const Scalar4 thermal_vel_mass = d_alt_vel[idx];
     vec3<Scalar> thermal_vel(thermal_vel_mass);
 
diff --git a/hoomd/mpcd/pytest/test_collide.py b/hoomd/mpcd/pytest/test_collide.py
index 83650811a8..aeaaae35a6 100644
--- a/hoomd/mpcd/pytest/test_collide.py
+++ b/hoomd/mpcd/pytest/test_collide.py
@@ -533,3 +533,39 @@ def test_rigid_mass_errors(
     # run simulation
     with pytest.raises(RuntimeError):
         sim.run(1)
+
+
+@pytest.mark.serial
+def test_rigid_nonparticipatory_zero_mass(small_snap, simulation_factory):
+    # create simulation
+    initial_snap = small_snap
+    if initial_snap.communicator.rank == 0:
+        initial_snap.particles.types = ["A", "B", "C"]
+        initial_snap.particles.mass[:] = [2]
+    sim = simulation_factory(initial_snap)
+    sim.seed = 5
+
+    rigid = hoomd.md.constrain.Rigid()
+    rigid.body["A"] = {
+        "constituent_types": ["B", "B", "C"],
+        "positions": [[-2.4, 0, 0], [2.4, 0, 0], [0, 1, 0]],
+        "orientations": [[1, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0]],
+    }
+    rigid.create_bodies(sim.state, masses={"A": np.array([1, 1, 0])})
+
+    sim.operations.integrator = hoomd.mpcd.Integrator(dt=0, rigid=rigid)
+    sim.operations.integrator.collision_method = (
+        hoomd.mpcd.collide.StochasticRotationDynamics(
+            period=1,
+            embedded_particles=hoomd.filter.Type("B"),
+            angle=90,
+            kT=1,
+        )
+    )
+
+    # run simulation
+    # velocities should all stay zero because there is no solvent to collide with
+    sim.run(1)
+    new_snap = sim.state.get_snapshot()
+    if new_snap.communicator.rank == 0:
+        assert np.array_equal(np.zeros((4, 3)), new_snap.particles.velocity)