UUID manager needs to be designed

[jngrad]

Shape-based constraints have a particle representation Particle ShapeBasedConstraint::part_rep which is default-constructed, therefore its particle id is always -1. This causes an issue when a real particle interacts with two shape-based constraints via the DPD interaction, which requires random numbers. Philox maps a tuple (counter, seed, salt, pid1, pid2) to a unique 4-tuple of 64bit integers. If a particle with id pid1 interacts with e.g. two parallel walls via the DPD interaction, pid2 will be -1 for both walls and the noise will be correlated.

Currently we work around this issue by incrementing the DPD thermostat counter. This has several issues:

• simulation reproducibility: the order by which the constraints are traversed in the short-range loop matters
• running system.integrator.run(0, recalc_forces=True) to update the forces alters the system state by incrementing the RNG counter
• we cannot use a single, monotonically increasing global counter

A couple ideas I and @KaiSzuttor came up with to solve the issue:

1. Instantiate "pseudo particles" such as the one used as a representation of shape-based constraints with the help of a function int get_maximal_pseudo_particle_id() that would work like the existing int get_maximal_particle_id() function. This function would be used to decrement the particle id, e.g. create a series -2, -3, -4, ..., because we must generate particle ids that live in the 32bit range while not colliding with particle ids from "real" particles (otherwise the correlation is not solved but moved to a different place). If such a function cannot be created, a global counter for "pseudo particles" could be used. Things to keep in mind:
   • nothing prevents the user from creating particles with values close to 2^32
2. Let the user instantiate real particles and use them as a parameter to objects that need a particle representation. For example, a wall constraint would be created with:

   p = system.part.add(pos=[0, 0, 0], type=2)
   shape_constraint = espressomd.constraints.ShapeBasedConstraint(shape=my_shape)
   constraint = system.constraints.add(shape=shape_constraint, particle_representation=p)

   instead of:

   shape_constraint = espressomd.constraints.ShapeBasedConstraint(shape=my_shape)
   constraint = system.constraints.add(shape=shape_constraint, particle_type=2)

   The constraint object would contain a shared pointer to the particle. Things to keep in mind:
   • the pseudo particle now lives in the cell system
   • particle slices will contain this particle, e.g. system.part[:].f now contains the constraint's particle "force" (it's unclear what its value should be)
   • integrators should not calculate its force, nor propagate it, nor use it to halt integration (steepest descent)
   • no particle coupling in LB
   • most members of the pseudo particle are meaningless (e.g. mass), possible issues with observables?

[jngrad]

Variant of option 1.:

3. Create a new salt DPDShapeBasedConstraint, so that particle ids of shape-based constraints can be incremented in the 32bit range using int get_maximal_pseudo_particle_id(). There is no collisions as long as the first pid in the key is always the shape-base constraint, for example.

In any case, we would have to expose the particle representation to the cython class to be able to checkpoint the particle id.

[KaiSzuttor]

In general, I'm not sure how big the issue is at all. You get correlated noise for a single particle and the constraint interaction.

Variant of option 1.:
Create a new salt DPDShapeBasedConstraint, so that particle ids of shape-based constraints can be incremented in the 32bit range using int get_maximal_pseudo_particle_id(). There is no collisions as long as the first pid in the key is always the shape-base constraint, for example.
In any case, we would have to expose the particle representation to the cython class to be able to checkpoint the particle id.

The meta problem is that we have no abstraction for the random number generation between particles and constraints.
The design of the RNG generation has to be rethought. We should not add hacks to circumvent this single issue here, otherwise we will hit the very same problem the next time we have a new interaction species.

[jngrad]

From the offline core team meeting:

4. Change RNG function to take a pair of uuids instead of a pair of particle ids. Add a new field in the Particle struct that stores a uuid. This uuid can be generated from a global counter. A prototype needs to be written in Python.