Apply a non-constant custom force to a specific particle

[reyhoon1998]

Hello everyone,

I'm trying to simulate a rod consisting of 4 particles that moves in the direction of its first particle.
So I generate a random orientation at initialization, and I use a custom force in HOOMD like this:

`def angle(pos_state):
pos = pos_state.particles.position
Lx, Ly = pos_state.configuration.box[0], pos_state.configuration.box[1]
# Unwrapped coordinates for all particles
x = pos[:, 0]
y = pos[:, 1]

theta = np.zeros(pos_state.particles.N) 
for i in head_tags:
    j = i + 1  
    dx_wr = x[i] - x[j]
    dy_wr = y[i] - y[j]

    if dx_wr > (Lx/2):
        dx_true = dx_wr - Lx
    elif dx_wr < - (Lx/2):
        dx_true = dx_wr + Lx
    else:
        dx_true = dx_wr

    if dy_wr > (Ly/2):
        dy_true = dy_wr - Ly
    elif dy_wr < - (Ly/2):
        dy_true = dy_wr + Ly
    else:
        dy_true = dy_wr

    theta[i] = np.arctan2(dy_true, dx_true)
return theta

Custom force for movement

class per_body_force(hoomd.md.force.Custom):
def init(self, CONSTANT):
super().init()
self.F = CONSTANT

def set_forces(self, timestep):
    snap = self._state.get_snapshot()
    if snap.communicator.rank != 0:
        return

    with self.cpu_local_force_arrays as arrays:
        theta = angle(snap)
        for i in head_tags:
            arrays.force[i] = self.F * np.array([np.cos(theta[i]), np.sin(theta[i]), 0.0])

custom_force = per_body_force(1)

LANGEVIN

langevin = hoomd.md.methods.Langevin(
filter=hoomd.filter.All(), kT=0.0, default_gamma=1
)

simulation.operations.integrator = hoomd.md.Integrator(
dt=0.01, methods=[langevin], forces=[bonding, rigidity, lj, custom_force]
)
gsd_writer = hoomd.write.GSD(
filename="move2body.gsd", trigger=hoomd.trigger.Periodic(1), mode="wb"
)

simulation.operations.writers.append(gsd_writer)

try:
simulation.run(100000)
finally:
gsd_writer.flush()
`

It doesn’t work as expected; I can see some unreasonable rotation, even though there is no randomness in the system.
I tried several ideas but couldn’t solve the issue.

1. Overall, can I use a custom force like this at all?

2.I suspect that HOOMD internally reorders particles during the simulation (for domain decomposition or optimization).
If that’s true, then my function, which depends on particle indices, may not refer to the same particle over time.
if so, is there a better way to implement this kind of directional force in HOOMD?

3. Another strange thing: I defined an unrelated custom action (triggered every step) that takes a snapshot and modifies something in it.When I enable that class, the problem with the custom force disappears! I’m really surprised; how could these two be related?

[joaander]

Overall, can I use a custom force like this at all?

Yes, you welcome to write whatever custom code you like.

I suspect that HOOMD internally reorders particles during the simulation...

Yes, the documentation makes a clear distinction between index order and tag order. The reordering increases performance by improving cache coherency and is employed even in serial execution.

• https://hoomd-blue.readthedocs.io/en/v5.4.0/tutorial/03-Parallel-Simulations-With-MPI/03-Accessing-System-Configurations-With-MPI.html#Modifying-Particle-Properties-With-Local-Snapshots
• https://hoomd-blue.readthedocs.io/en/v5.4.0/hoomd/data/particlelocalaccessbase.html#hoomd.data.ParticleLocalAccessBase.tag

When I enable that class, the problem with the custom force disappears! I’m really surprised; how could these two be related?

I do not know.

I think you are over-complicating the solution to your stated goal. To drive an active rod, you can use md.force.Active to apply a force on the first particle. You provide the active force in the local coordinate system of that particle (e.g. (0,0,1)) and force.Active will rotate that force into the reference frame of that particle and apply it on each timestep.

Here is a work-in-progress tutorial on active matter simulations that explains more: https://github.com/glotzerlab/hoomd-examples/tree/activity-tutorial/11-Active-Particles-In-MD

[reyhoon1998]

Thank you so much for your help. I read all documentation you sent, and can't wait to see the complete tutorial.