Particles without user-defined type are not tracked by reactions

[jngrad]

When inserting particles of type 0 with the following syntax:

system.part.add(pos=[1, 2, 3])

the particle is created with a default-constructed value for the type, in this case 0. Since the set_particle_type() callback is not called, the particle id is not added to the global type map used by reaction methods, and therefore won't be involved in any reaction where one of the reactants has type 0. This affects all ESPResSo versions since 4.0.0.

MWE for 4.2.0 and 4.3-dev:

import espressomd.reaction_methods
import numpy as np
np.random.seed(seed=42)
system = espressomd.System(box_l=3 * [35.0], time_step=0.02)
system.cell_system.skin = 0.4
# 2A + 3B <-> 4C
type_A, type_B, type_C = (0, 1, 2)
nu_A, nu_B, nu_C = (-2, -3, 4)
nu_bar = sum([nu_A, nu_B, nu_C])
K = 1000.
RE = espressomd.reaction_methods.ReactionEnsemble(
    kT=1., exclusion_range=0., seed=4)
RE.set_non_interacting_type(type=3)
RE.add_reaction(
    gamma=K * (0.027)**nu_bar, reactant_types=[type_A, type_B],
    reactant_coefficients=[abs(nu_A), abs(nu_B)],
    product_types=[type_C],
    product_coefficients=[nu_C],
    default_charges={type_A: 0, type_B: 0, type_C: 0})
N0 = 50
for i in range(N0):
    for j in range(abs(nu_B)):
        system.part.add(pos=np.random.random(3) * system.box_l, type=type_B)
    for j in range(abs(nu_A)):
        system.part.add(pos=np.random.random(3) * system.box_l) # missing type=type_A
RE.reaction(reaction_steps=2000)
print(f"n_A = {system.number_of_particles(type=type_A)}")
print(f"n_B = {system.number_of_particles(type=type_B)}")
print(f"n_C = {system.number_of_particles(type=type_C)}")

MWE for 4.1.4:

import espressomd.reaction_ensemble
import numpy as np
np.random.seed(seed=42)
system = espressomd.System(box_l=3 * [35.0], time_step=0.02)
system.cell_system.skin = 0.4
# 2A + 3B <-> 4C
type_A, type_B, type_C = (0, 1, 2)
nu_A, nu_B, nu_C = (-2, -3, 4)
nu_bar = sum([nu_A, nu_B, nu_C])
K = 1000.
RE = espressomd.reaction_ensemble.ReactionEnsemble(
    temperature=1., exclusion_radius=0., seed=4)
RE.set_non_interacting_type(3)
RE.add_reaction(
    gamma=K * (0.027)**nu_bar, reactant_types=[type_A, type_B],
    reactant_coefficients=[abs(nu_A), abs(nu_B)],
    product_types=[type_C],
    product_coefficients=[nu_C],
    default_charges={type_A: 0, type_B: 0, type_C: 0})
N0 = 50
for i in range(N0):
    for j in range(abs(nu_B)):
        system.part.add(pos=np.random.random(3) * system.box_l, type=type_B)
    for j in range(abs(nu_A)):
        system.part.add(pos=np.random.random(3) * system.box_l) # missing type=type_A
RE.reaction(reaction_steps=2000)
print(f"n_A = {system.number_of_particles(type=type_A)}")
print(f"n_B = {system.number_of_particles(type=type_B)}")
print(f"n_C = {system.number_of_particles(type=type_C)}")

Output:

n_A = 0
n_B = 150
n_C = 0

Expected output:

n_A = 40
n_B = 60
n_C = 120

Possible patch for 4.3-dev:

diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp
index 70b9330a93..19e3f15323 100644
--- a/src/script_interface/particle_data/ParticleHandle.cpp
+++ b/src/script_interface/particle_data/ParticleHandle.cpp
@@ -484,4 +484,7 @@ void ParticleHandle::do_construct(VariantMap const &params) {
       }
     }
+    if (params.count("type") == 0) {
+      do_set_parameter("type", 0);
+    }
   }
 }

Possible patch for 4.2.0:

diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx
index 3dc80634b..964143eee 100644
--- a/src/python/espressomd/particle_data.pyx
+++ b/src/python/espressomd/particle_data.pyx
@@ -1814,4 +1814,7 @@ Set quat and scalar dipole moment (dipm) instead.")
         pid = p_dict.pop("id")
 
+        if "type" not in p_dict:
+            p_dict["type"] = 0
+
         if p_dict != {}:
             self.by_id(pid).update(p_dict)

Possible patch for 4.1.4:

diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx
index 4e7799a29..c93b752fe 100644
--- a/src/python/espressomd/particle_data.pyx
+++ b/src/python/espressomd/particle_data.pyx
@@ -1912,4 +1912,7 @@ Set quat and scalar dipole moment (dipm) instead.")
         del P["id"]
 
+        if "type" not in P:
+            P["type"] = 0
+
         if P != {}:
             self[id].update(P)