Active force problem.

[Ye Zhang]

HOOMD-V3

I'm using `md.force.Active and md.force.custom_force to simulate active matter system with brownian integrater

md.force.active_force

active = hoomd.md.force.Active(filter=hoomd.filter.All())
active.active_force['A'] = (10,0,0)
integrator.forces.append(active)
rotational_diffusion_updater = active.create_diffusion_updater(rotational_diffusion=1,trigger=1)
sim.operations+=rotational_diffusion_updater

I want to change the active force during the simulation, so according to the document on active force,

https://hoomd-blue.readthedocs.io/en/stable/module-md-update.html#hoomd.md.update.ActiveRotationalDiffusion

md.force.custom_force

​
class CustomActiveForce_2D(hoomd.md.force.Custom):
 def __init__(self,f_array,theta_array,rotation_diff):
   super().__init__(aniso=True)
   self.f_array = f_array.reshape((-1,1))
   self.theta_array = theta_array.reshape((-1,1))
   self.rotation_diff = rotation_diff
   
 def update_force(self):
   self.theta_array+=np.random.normal(0,1,self.theta_array.shape)*(
     2*self.rotation_diff*integrator.dt)**0.5
   self.active_f = self.f_array*np.column_stack([np.cos(self.theta_array),
                   np.sin(self.theta_array),
                   np.zeros(self.theta_array.shape)])
   
 def set_forces(self,timestep):
   self.update_force()
   with self.cpu_local_force_arrays as arrays:
     arrays.force[:] =self.active_f
     pass
   
active_N = snapshot.particles.N
f_array =  10*np.ones((active_N,1))
custom_active =CustomActiveForce_2D(f_array = f_array,
                        theta_array= 2*np.pi*(np.random.random((snapshot.particles.N,1))-0.5),
                        rotation_diff=1)
integrator.forces.append(custom_active)
integrator.integrate_rotational_dof = True

But the result of the two methods is different when keeping the constants the same.

It's there anything wrong with my code?

Thanks for your help~

[Joshua Anderson]

Ye Zhang,

I'm not surprised the two implementations provide different results. You are generating different random numbers than HOOMD's. I can't say whether anything is wrong with your code. Only you can make such a judgement.

Here is HOOMD's implementation if you want to compare:
https://github.com/glotzerlab/hoomd-blue/blob/8fc5b26a4a28d42c86afc3740437aea66578af19/hoomd/md/ActiveForceCompute.cc#L249-L275
------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/8d0dbab8-f25f-40bd-a92a-1c2d20c4bd8fn%40googlegroups.com.
> <activeforce_question_HOOMD-V3.md>

[Ye Zhang]

thanks for your help~ I'm trying to solve this problem

[Joshua Anderson]

Problem descriptions must describe the behavior you get vs the expected behavior. Without such a description, we can only guess at what you intend to accomplish with your code and only guess at what your code actually does. Normally, people on the mailing list don't waste time offering such guesses.

I did think of another guess to offer: You store theta in index order. HOOMD reorders particles in memory, so theta[i] will apply to particle with tag t1 on one timestep and tag t2 on the next. To completely replicate hoomd.md.force.Active, you should use the particle orientation to track the direction of the force vector. This enables support for MPI domain decomposition as well.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/6bf827ea-3ee4-4309-ac8d-2e29d7668969n%40googlegroups.com.

[Ye Zhang]

thanks for your help~

[Ramalingam Kailasham]

Hi All

I recently got started with using HOOMD to simulate active Brownian particles. I agree with Joshua's recommendation that using the particle orientation (quaternion) to update the active force is the correct way to go about implementing active forces in a custom manner. 

With regards to Ye Zhang's code, the formulae appear correct. However, the custom force is sampling "theta" from a random distribution at each timestep, while one would want the updated theta to remember the value of the theta from the previous timestep. This is the reason that working with the particle orientations is preferred. 

Hope that helps

Regards