Add support for writing rigid bodies when converting to HOOMD formats.

[chrisjonesBSU]

rigid IDs are ints assigned to particles/atoms that signal which rigid body they belong to which is needed when recording rigid body information in HOOMD Frames and/or GSD files. The mBuild hoomd writers handled this, but we still need to implement them here.

In mBuild, rigid_id is a property at the particle level, so we can treat it the same here, at the site level.

This is a draft for now. Remaining TODOs are:

• handle cases where some molecules are rigid and others aren't
• throw an error if rigid particle type is found in the topology already

[CalCraven]

A remaining question to figure out: What if only some molecules are rigid in a hoomd simulation? How does the list of rigid body indices change?
https://hoomd-blue.readthedocs.io/en/v3.6.0/tutorial/06-Modelling-Rigid-Bodies/00-index.html

[CalCraven]

Was playing around with the rigid bodies in HOOMD. Got a working simulation with some rigid particles and some non-rigid particles. Looks like the body argument, shown now in line, get's a tag of -1 if the particle is not rigid. The rigid particles get the rigid ID for that particle, and the rigid COM particle, named "dimer" here, also get's that same ID.

[chrisjonesBSU]

The simulations with rigid bodies seem to be working as expected. Here is a Gist of the notebook I'm using to run them if you want to try it.

A couple thoughts.

• to_gsd_snapshot() and to_hoomd_snapshot() now return a dictionary of rigid body info (along with the snapshot and reference values dictionary) that has information needed to create the hoomd.md.constrain.Rigid object. After thinking about it more, we could just create hoomd.md.constrain.Rigid, populate the dictionary and return it.

This technically breaks the old behavior, since these methods now return 3 objects instead of 2. Also, being required to create a variable or placeholder for this in workflows that don't use rigid bodies might be awkward. We could add a bool parameter to these methods. Something like create_hoomd_constraint that if toggled True will also create and return the hoomd.md.constrain.Rigid object. Or, we could do this automatically, and the return statement used depends on if there are rigid bodies or not. I'm not sure which is preferred or best practice. Ultimately, I think it would be very helpful if the hoomd.md.constrain.Rigid can be made at the same time as the snapshot.

• Should to_hoomd_forcefield() have some parameters/checks for rigid bodies as well?

the non-bonded force objects need to include parameters between the rigid particles and their constituent particles. Also, the neighbor list exclusions are different for rigid body simulations than those without. The gist handles these manually after calling to_hoomd_forcefield(). Maybe this method should make these changes if the topology contains rigid bodies.

[chrisjonesBSU]

I've been testing this out a bit. It works with multiple types of rigid molecules (e.g., a system with rigid benzene and rigid ethane), and it works with a mix of rigid and non-rigid molecules (e.g., rigid benzene and flexible ethane), BUT the rigid molecules have to come first in the compound/topology hierarchy.

We can discuss more on a dev call, but there are still some questions.

I think the gsd writers should also build the hoomd.md.constrain.Rigid object that is needed to run rigid simulations, since all of the information needed is right there when creating the snapshot, but:

Do we make it so the gsd writers always return 3 things now (breaking change) or have the objects returned depend on whether or not rigid bodies exist (not necessarily a breaking change)

So:

return snapshot, ref_values, rigid_constrain

where rigid_constrain is just None if there aren't rigid bodies

or

if rigid_bodies:
    return snapshot, ref_values, rigid_constrain
else:
    return snapshot, ref_values

The other thing to think about is what behavior we want in to_hoomd_forcefield. If the user is using rigid bodies do we:

• Make the necessary changes to the neighbor list (update exclusions to include "body")
• Automatically handle the LJ params between rigid particles and constituent atoms (i.e., add parameters needed to ignore these interactions)
• For systems where every molecule is rigid, do we remove any intramolecular interactions (bonds, angles, dihedrals, impropers) from the forces returned, or just skip their parsing methods to begin with?

Or, do we just return the complete forcefield as is, and leave it to the user to handle these things? In this case, we could make sure to create a thorough tutorial in the documentation, and provide some examples in a notebook.

One idea is possibly having a wrapper around the hoomd writers specific to rigid body stuff, so that the existing writers stay true to form as much as possible. Maybe something like from external.convert_hoomd import create_rigid_body_sim that calls to_gsd_snapshot and to_hoomd_forcefield under the hood, and makes the above changes.