Active Forces on Rigid Bodies

Overview

Questions

• How do we correctly model active forces on an arbitrary rigid body in HOOMD-blue?

Objectives

• Define a rigid body as a central particle and a set of constituent particles in a lattice.

• Understand how to model driven active matter systems by applying an active force to the rigid body.

Boilerplate Code

import io
import itertools

import gsd.hoomd
import hoomd
import numpy
import PIL

The render_rigid_body function in the next (hidden) cell will render the rigid body using fresnel. Find the source in the hoomd-examples repository.

Applying active forces to rigid bodies

This tutorial demonstrates how to simulate active Brownian particles made up of rigid bodies using HOOMD-blue. Accurately applying active forces to rigid bodies is essential for realistic physical simulations. Active forces and torques can arise from self-propulsion or self-rotational mechanisms, external fields, or internal actuation, and they strongly influence collective dynamics and emergent behavior.

In the previous tutorials, we applied forces on the central particle of rigid bodies. Here we do the same thing, but with a linear active force to achieve self-propulsion of our rigid body. Alternatively, adding active torque to the central particle of our rigid body would simulate self-rotation.

One could alternatively apply an active force on a constituent particle that is not located on the center of mass. Any active torque applied not in line with the radius will induce an active torque. It is recommended to avoid doing so. To achieve this effect, one should instead consider applying an active torque on your system as below unless there is a good reason for not doing so.

For a full tutorial on modeling rigid bodies please see the Modeling Rigid Bodies tutorial.

Here we give an example of a system where only active torques are applied. Some relevant example systems include the development of starfish embryos or using a bacterial bath to assemble active colloids.

To model this system, we will start by seting up a similar simulation to the previous tutorials.

Set up the Simulation

We will first initialize a system of rigid bodies with a given particle density in an array. We create an initial configuration of dimers with a frame containing both particle types.

m = 4
N_particles = m**3
spacing = 5
K = numpy.ceil(N_particles ** (1 / 3))
L = K * spacing
x = numpy.linspace(-L / 2, L / 2, int(K), endpoint=False)
position = list(itertools.product(x, repeat=3))
position = numpy.array(position) + [spacing / 2, spacing / 2, spacing / 2]

frame = gsd.hoomd.Frame()
frame.particles.N = N_particles
frame.particles.position = position[0:N_particles, :]
frame.particles.typeid = [0] * N_particles
frame.configuration.box = [L, L, L, 0, 0, 0]
frame.particles.types = ["dimer", "A"]
frame.particles.mass = [2] * N_particles
frame.particles.orientation = [(1, 0, 0, 0)] * N_particles

Compute the moment of inertia of the dimer. Computing this is essential to accurate rotational updating of the rigid body.

particle_mass = 1
particle_radius = 1
dimer_positions = [[1.2, 0, 0], [-1.2, 0, 0]]
I = numpy.zeros(shape=(3, 3))
for r in dimer_positions:
    I += particle_mass * (numpy.dot(r, r) * numpy.identity(3) - numpy.outer(r, r))
frame.particles.moment_inertia = [I[0, 0], I[1, 1], I[2, 2]] * N_particles

Write the rigid centers to a GSD file for later use:

with gsd.hoomd.open(name="lattice.gsd", mode="w") as f:
    f.append(frame)

Configure the simulation to run on the CPU using this snapshot, as done in previous tutorials. Here, we also define the constituent particle properties in body coordinates for the rigid body type. After creating the rigid bodies, visualize the resulting configuration:

simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed=1)
state = simulation.create_state_from_gsd(filename="lattice.gsd")

rigid = hoomd.md.constrain.Rigid()
rigid.body["dimer"] = {
    "constituent_types": ["A", "A"],
    "positions": dimer_positions,
    "orientations": [(1.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0)],
}
rigid.create_bodies(simulation.state)
render(simulation.state.get_snapshot())

Apply the OverdampedViscous integrator and set kT=0 to make active translation more obvious. We want to allow our rigid bodies to rotate, so it’s important to integrate the rotational degrees of freedom. Thus, we set integrate_rotational_dof = True unlike in the previous tutorial notebook.

Importantly, note that when defining the integrator, we have added the argument rigid=rigid. This is the argument that tells the HOOMD integrator to treat the particles defined by the Rigid object as a single rigid body, so it integrates only the center particle’s translation and rotation and then updates the constituent particles’ positions/orientations to maintain the fixed internal geometry.

filter_all = hoomd.filter.All()
filter_free = hoomd.filter.Rigid(("center", "free"))
filter_dimer = hoomd.filter.Type(["dimer"])

gamma = 1

overdamped = hoomd.md.methods.OverdampedViscous(filter=filter_free)
overdamped.gamma.default = gamma
overdamped.gamma_r.default = [gamma / 3, gamma / 3, gamma / 3]


integrator = hoomd.md.Integrator(dt=0.001, rigid=rigid, integrate_rotational_dof=True)
simulation.operations.integrator = integrator
integrator.methods.append(overdamped)

Apply our active force

When applying an active force to a rigid body, it matters which particle the active force is applied to. This can be applied to the central particle of the rigid body (dimer in the present case) or the body’s constituent particles. Regardless of which particles you choose, the central particle holds all the information provided regarding the force. Only one active force should be applied to the central particle of the rigid body. In this case this particle is dimer. We recommend that this central particle is the center of mass of the rigid body. Not doing so is possible, but it is not particularly recommended.

Notice that we do not create a rotational_diffusion_updater here unlike in previous tutorials of active forces. We do so only for point particles. This is because, in HOOMD-blue, a rigid body’s orientation is one of the system’s dynamical degrees of freedom, and the integrator updates it by solving the rotational equations of motion. If we were to apply a rotational_diffusion_updater to a rigid body, we would be manually changing its orientation in addition to the integrator’s rotational dynamics. This would effectively add an extra source of rotational noise that is not part of the rigid-body equations of motion and would generally be unphysical.

Below we apply an active force in the x and y body coordinates solely because it provides easier visualization for the sake of this tutorial.

activity = 0.3  # x and y force component

active = hoomd.md.force.Active(filter=filter_dimer)
active.active_force["dimer"] = (activity, activity, 0)
active.active_torque["dimer"] = (0, 0, 0)

integrator.forces.append(active)

Write simulation data to a GSD trajectory file every 100 steps.

logger = hoomd.logging.Logger(["particle", "constraint"])
gsd_oper = hoomd.write.GSD(
    trigger=hoomd.trigger.Periodic(100),
    filename="active.gsd",
    logger=logger,
    mode="wb",
    dynamic=["momentum", "property"],
    filter=filter_all,
)

simulation.operations += gsd_oper

Visualize the system’s initial configuration to compare to later trajectory frames.

simulation.run(0)
render(simulation.state.get_snapshot())

Run the simulation. Since we have added no thermal fluctuations to this simulation, any observed particle motion results from the active propulsion force applied to the rigid bodies. Here, the entire rigid body is propelled by the force applied to its center of mass because the active_force is applied to the central body particle dimer. The body’s orientation and position evolves over time, and the constituent particles follow accordingly. Although the force is explicitly applied only to the central dimer particle, the entire rigid body translates collectively, as expected. This confirms that the active force is correctly coupled to the rigid body and propagated to its constituent particles.

simulation.run(5000)
gsd_oper.flush()
render(simulation.state.get_snapshot())

Apply an Active Torque

Finally, we will demonstrate the application of an active torque on this system. This again requires that we set the argument in the simulation integrator integrate_rotational_dof=True. For ease of visualization of the active torque, we will first remove the linear active force so that the center of mass of each dimer will remain stationary.

integrator.forces.remove(active)

We now apply an active torque analogously to the linear active force.

tau = 0.3

active = hoomd.md.force.Active(filter=filter_dimer)
active.active_force["dimer"] = (0, 0, 0)
active.active_torque["dimer"] = (0, 0, tau)

integrator.forces.append(active)

simulation.run(1000)
gsd_oper.flush()
render(simulation.state.get_snapshot())

All observed change in orientation is due to the applied active torque on our dimers since we have added no thermal noise to this system. We can verify by visual inspection that each of our dimers are indeed rotated by an active torque as individual rigid bodies.