Support rigid bodies in MPCD

[mzbush]

Description

Rigid bodies are currently incompatible with MPCD collision methods because:

1. The collision method needs to read the velocities of constituent particles, but they are not currently set.
2. The collision method updates the velocity of constituent particles in a way that is inconsistent with rigid motion.

We should add support for rigid bodies to the collision methods.

Proposed solution

1. ForceComposite should set the velocity of constituent particles based on the linear and angular velocity of the rigid body when it sets their positions and orientations.

   • HOOMD's angular momentum quaternion can be used to compute the angular velocity in the body frame.

   • The body's angular velocity can be used to get constituent particle velocities in the body frame relative to the center-of-mass translation.

   • These velocities can then be rotated to the space frame and have the center-of-mass translational velocity added.

   • Question: should this always be done, or should it be opt in?

2. The CollisionMethod should store the net change in space-frame linear and angular momentum for each rigid body as a result of collisions on constituent particles.

   • The changes can be computed by the base class taking the difference in velocity of constituent particles before/after collision or by recording the change during each collision. The first option is more generic, but the second option is probably faster because it has less read/write transactions.

   • In MPI simulations, we will need to reduce these contributions onto the rank that owns the central particle, then broadcast the result.

   • These changes will be used to update the velocity and angular momentum of the central particle for the body, then reset the constituent velocities.

Additional context

No response

[joaander]

Point 1 is briefly discussed in #278. Note that implementing it with MPI will require communicating the ghost angular momentum so that all ghost constituents have the needed information to compute their velocity locally. For consistency with the rest of HOOMD-blue, the constituent velocities should always be set. I think there is a flags mechanism to avoid communicating angular momentum when it is not needed.

[mphoward]

This is a good point: DPD (#278) has a similar issue as MPCD, in that the momentum transferred to constituent particles needs to be passed to the central particle to enforce the rigid constraint. We should consider this when we design our solution for MPCD so that maybe the same code can also be used for DPD.

Our plan is to set the velocities of the constituent particles first (making sure it works CPU, GPU & MPI), then we'll figure out how to address the momentum transfer part.

[joaander]

DPD applies forces to the constituents which will then be transferred as a net force/torque to the central particle by the normal means. The problem with using DPD and rigid bodies is that the constituents do not have a velocity set and the DPD thermostat depends on the relative velocity of the particle pair. Unless I misunderstood your solution 1 would solve this by setting the constituent velocities.

[mphoward]

DPD applies forces to the constituents which will then be transferred as a net force/torque to the central particle by the normal means. ... Unless I misunderstood your solution 1 would solve this by setting the constituent velocities.

No, you are right! I had in my head that there was a separate path taken for the thermostatting part, but since everything is done as a force, simply setting the velocities beforehand should be sufficient. We can then take the optimal path for MPCD to transfer the momentum, which is not done using a force.

[mphoward]

While working on this @mzbush realized that for MPCD we will need to be able to set the masses of the constituent particles, but they are currently ignored because the constituents don't have velocities.

One way to do this would be to add a masses key to the per-body properties of constituent particles. To preserve backward compatibility, this might need to default to something sensible.

I noticed, though, that charges of constituent particles are currently optionally set through create_bodies rather than the per-body properties. That would be another option that could work, but I was unsure why this was handled in this way rather than the mechanism above.

@joaander what do you advise?

[joaander]

Constituent particles are merely interaction sites for potentials. That is why they only have a position, orientation, and type defined in the local reference frame of the body. They do not have mass because they have no degrees of freedom (and therefore no inertia) on their own.

Due to the fact that HOOMD-blue stores constituent particles in the same array with the body centers, there is data storage for mass, velocity, angular momentum, charge, etc... With the exception of charge (and velocity as discussed above), none of the other fields are used in the computation of potentials. create_bodies accepts charge as a convenience for users wanting to place constituent particles in the state with the given charges (it was a commonly used feature in HOOMD-blue 2.x). However, create_bodies is completely optional, you can place the constituents manually in your state and give them whatever properties you want. position, orientation, and type will be overridden by constrain.Rigid, but all others should be left alone.

I realize this behavior is confusing. I don't have a solution to fix that in HOOMD-blue. I would prefer not to add mass to create_bodies. It will encourage the false sense that somehow the mass of the body is the sum of the masses of the constituents.

[mphoward]

Thanks! This makes sense for the current design. Now that we are setting constituent velocities, mass is a bit like type from an implementation perspective: it doesn't depend on local vs. global frame, but we need to set something as the w element of the velocity when we overwrite those for constituent particles (like how type gets written into the w element of updated position). I was wondering if it might make sense to then include masses as part of the body definition.

We are currently reading the value of the mass for the constituent particles, then writing it back with the updated velocity. We can continue to do that, and it will mean that mass needs to be manually set if the user wants to use it. We can experiment with how that works for us in practice and reconsider if it is too clunky.

[joaander]

Yes, because velocity and mass are linked you will need to read the current mass in and write it back out to maintain the current data model. I agree that it will be clunky. Perhaps the rigid body tutorials can be improved to show more clearly how users should initialize parameters like mass for their constituents. I find that users often don't understand what create_bodies does for them and use it in a variety of unintended ways. Most help sessions with users trying to use rigid bodies end when we find either a missing call to create_bodies or extra unneeded calls to create_bodies that produce unexpected behavior.

[mphoward]

Perhaps the rigid body tutorials can be improved to show more clearly how users should initialize parameters like mass for their constituents.

We will definitely need to do this for our MPCD tutorial on rigid bodies so it wouldn't be hard to add similar general instructions to the main rigid body tutorial as well. Mass and charge would both be good examples.