Active torque makes FENEWCA bond incorrect

[kizzhang]

Description

I encountered this problem as I simulated active chiral particles chained into a polymer using FENEWCA potential. The simulation does not produce any FENE out-of-bound error, but apparently the bonds are way longer than they should be (10 vs 0.97).

Script

import hoomd
import gsd.hoomd

############################# GET PARAMETERS #################################

para_arr = sys.argv
seed_num = 1213
Pe_arg = 100.0
print('Pe_arg is ',Pe_arg)    
N_particles = 200

######################## Read Ring Polymer Configuration #########################
traj = gsd.hoomd.open('ring.gsd')
snapshot = traj[-1]

########################### RUN THE SIMULATION ###############################
cpu = hoomd.device.CPU()
sim = hoomd.Simulation(device=cpu, seed=seed_num)
sim.create_state_from_snapshot(snapshot)
print('Seed number is ',seed_num)

###########################  CREATE FORCE   ##################################

if Pe_arg <=1:
    kbt = 1
else:
    kbt = 1/ Pe_arg
print("kbt is = " + str(kbt))

f_sigma = 1.0
f_epsilon = 1.0
R_0 = 1.5*f_sigma

k_fene = 60 * f_epsilon/f_sigma**2

fenewca = hoomd.md.bond.FENEWCA()
fenewca.params['A-A'] = dict(k= k_fene, r0 = R_0, epsilon= f_epsilon, sigma= f_sigma,
                             delta=0.0)

nl = hoomd.md.nlist.Tree(buffer = 0.4)

lj = hoomd.md.pair.LJ(nl)
lj.params[('A', 'A')] = dict(epsilon=f_epsilon, sigma=f_sigma)
lj.r_cut[('A', 'A')] = 2**(1/6)*f_sigma

######################### CREATE DIFFUSION  ################################

gamma = np.sqrt(3)
dis_diff = kbt/gamma 
f_sigma = 1.0
rot_diff = 3*dis_diff/f_sigma**2
print("Rotational Diffusion is = ",rot_diff)

F_active = 1/np.sqrt(3)
print('F_active = ', F_active)
    
Pe = np.sqrt(F_active**2 * 3)*f_sigma/kbt
print("Pe is = "+str(Pe) +"; Pe_arg is = "+ str(Pe_arg))

brownian = hoomd.md.methods.Brownian(filter=hoomd.filter.All(),kT = kbt)
brownian.gamma.default = gamma
brownian.gamma_r.default = [0, 0, 0]

active = hoomd.md.force.Active(filter=hoomd.filter.All())
active.active_force['A'] = (F_active,F_active,F_active)
active.active_torque['A'] = (1,0,0)

integrator = hoomd.md.Integrator(dt=1e-5, methods=[brownian], forces=[fenewca,lj,active])#,integrate_rotational_dof=True)
rotational_diffusion_updater = active.create_diffusion_updater(trigger=hoomd.trigger.After(0), rotational_diffusion = rot_diff)
sim.operations += rotational_diffusion_updater

sim.operations.integrator = integrator

sim.run(1e5)
gsd_writer = hoomd.write.GSD(filename='out.gsd',
                             mode='xb',
                             trigger = hoomd.trigger.Periodic(int(5e4)),
                            )
sim.operations.writers.append(gsd_writer)
integrator.dt = 1e-4
sim.run(1e7)
gsd_writer.flush()

Input files

ring.txt

Output

As illustrated in the figure, the bond is over-stretched while FENEWCA calculation seems overridden. This happens with torques in Constant.Force as well. Things are normal if not including torque.

Expected output

Expecting the bond to be restricted by FENEWCA potential or output bond out-of-bound error.

Platform

Linux, CPU

Installation method

Conda-forge package

HOOMD-blue version

5.1.0

Python version

3.10.13

[joaander]

Thanks for the complete report. I'll try running your script when I have the time to. One quick question first:
Do you see the bonds break with an active torque set and the default integrate_rotational_dof=False? Or only when integrate_rotational_dof=True?

I don't know how the torque could have any effect when integrate_rotational_dof=False, as the torques will be set but have no effect on the motion of the particles.

Also, what is your intention to use both an active torque and a rotational diffusion updater? Both will change the orientations of the particles in different ways that do not interoperate. Normally, one would use a rotational diffusion updater with point particles and an active torque on anisotropic particles that drive themselves to always rotate in the same direction.

[kizzhang]

Thanks for the complete report. I'll try running your script when I have the time to. One quick question first: Do you see the bonds break with an active torque set and the default integrate_rotational_dof=False? Or only when integrate_rotational_dof=True?

I don't know how the torque could have any effect when integrate_rotational_dof=False, as the torques will be set but have no effect on the motion of the particles.

Also, what is your intention to use both an active torque and a rotational diffusion updater? Both will change the orientations of the particles in different ways that do not interoperate. Normally, one would use a rotational diffusion updater with point particles and an active torque on anisotropic particles that drive themselves to always rotate in the same direction.

I appreciate your timely response.

I ran the simulation with the default option ( integrate_rotational_dof=False).

I am trying to replicate this simulation (I noticed that Eqn. 2ab is quite different from representation in doc. I looked at the activeCompute.cc file, the rotated vector is represented as quaternion product, which is also different from the doc too.) I tried to include "chirality", where the paper used Gaussian noise with a non-zero mean (It is simple to implement chirality in 2D, i.e. $d\theta/dt=\omega+\xi$). I understand that what I am doing now is not the correct way to implement Eqn. 2ab from the paper. But this bug happened when I was playing around with HOOMD.

[joaander]

I attempted to run your script. The first error I get is:

  File "/Users/joaander/test/ring-active-torque/test.py", line 6, in <module>
    para_arr = sys.argv
               ^^^
NameError: name 'sys' is not defined. Did you forget to import 'sys'?

para_arr does not appear to be used so I deleted the line. The next problem is that ring.txt is a text file, not a .gsd file. Without ring.gsd, I cannot reproduce the behavior you describe.

I did confirm that Brownian will ignore all torque values when integrate_rotational_dof == False:

hoomd-blue/hoomd/md/TwoStepBD.cc

Lines 154 to 246 in 412c2ec

	if (m_aniso)
	{
	unsigned int type_r = __scalar_as_int(h_pos.data[j].w);
	Scalar3 gamma_r = h_gamma_r.data[type_r];
	if (gamma_r.x > 0 || gamma_r.y > 0 || gamma_r.z > 0)
	{
	vec3<Scalar> p_vec;
	quat<Scalar> q(h_orientation.data[j]);
	vec3<Scalar> t(h_torque.data[j]);
	vec3<Scalar> I(h_inertia.data[j]);
	bool x_zero, y_zero, z_zero;
	x_zero = (I.x == 0);
	y_zero = (I.y == 0);
	z_zero = (I.z == 0);
	Scalar3 sigma_r
	= make_scalar3(fast::sqrt(Scalar(2.0) * gamma_r.x * currentTemp / m_deltaT),
	fast::sqrt(Scalar(2.0) * gamma_r.y * currentTemp / m_deltaT),
	fast::sqrt(Scalar(2.0) * gamma_r.z * currentTemp / m_deltaT));
	if (m_noiseless_r)
	sigma_r = make_scalar3(0, 0, 0);
	// original Gaussian random torque
	// Gaussian random distribution is preferred in terms of preserving the exact math
	vec3<Scalar> bf_torque;
	bf_torque.x = NormalDistribution<Scalar>(sigma_r.x)(rng);
	bf_torque.y = NormalDistribution<Scalar>(sigma_r.y)(rng);
	bf_torque.z = NormalDistribution<Scalar>(sigma_r.z)(rng);
	if (x_zero)
	{
	bf_torque.x = 0;
	t.x = 0;
	}
	if (y_zero)
	{
	bf_torque.y = 0;
	t.y = 0;
	}
	if (z_zero)
	{
	bf_torque.z = 0;
	t.z = 0;
	}
	// use the damping by gamma_r and rotate back to lab frame
	// Notes For the Future: take special care when have anisotropic gamma_r
	// if aniso gamma_r, first rotate the torque into particle frame and divide the
	// different gamma_r and then rotate the "angular velocity" back to lab frame and
	// integrate
	bf_torque = rotate(q, bf_torque);
	if (D < 3)
	{
	bf_torque.x = 0;
	bf_torque.y = 0;
	t.x = 0;
	t.y = 0;
	}
	// do the integration for quaternion
	q += Scalar(0.5) * m_deltaT * ((t + bf_torque) / vec3<Scalar>(gamma_r)) * q;
	q = q * (Scalar(1.0) / slow::sqrt(norm2(q)));
	h_orientation.data[j] = quat_to_scalar4(q);
	if (m_noiseless_r)
	{
	p_vec.x = t.x / gamma_r.x;
	p_vec.y = t.y / gamma_r.y;
	p_vec.z = t.z / gamma_r.z;
	}
	else
	{
	// draw a new random ang_mom for particle j in body frame
	p_vec.x = NormalDistribution<Scalar>(fast::sqrt(currentTemp * I.x))(rng);
	p_vec.y = NormalDistribution<Scalar>(fast::sqrt(currentTemp * I.y))(rng);
	p_vec.z = NormalDistribution<Scalar>(fast::sqrt(currentTemp * I.z))(rng);
	}
	if (x_zero)
	p_vec.x = 0;
	if (y_zero)
	p_vec.y = 0;
	if (z_zero)
	p_vec.z = 0;
	// !! Note this isn't well-behaving in 2D,
	// !! because may have effective non-zero ang_mom in x,y
	// store ang_mom quaternion
	quat<Scalar> p = Scalar(2.0) * q * p_vec;
	h_angmom.data[j] = quat_to_scalar4(p);
	}

Brownian will also zero all torques applied to particle degrees of freedom where the moment of inertia is 0. Therefore, I still don't see how changing the active torque value in your script will change any behavior in the simulation. I will need a complete minimal working example to investigate further.

The active force and rotational diffusion documentation clearly state how they operate on the quaternions. From the active force docstring:

    `Active` computes an active force and torque on all particles selected by
    the filter:

    .. math::

        \vec{F}_i = \mathbf{q}_i \vec{f}_i \mathbf{q}_i^* \\
        \vec{\tau}_i = \mathbf{q}_i \vec{u}_i \mathbf{q}_i^*,

    where :math:`\vec{f}_i` is the active force in the local particle
    coordinate system (set by type `active_force`) and :math:`\vec{u}_i`
    is the active torque in the local particle coordinate system (set by type
    in `active_torque`.

And the rotational diffusion updater states:

    `ActiveRotationalDiffusion` works directly with `hoomd.md.force.Active` or
    `hoomd.md.force.ActiveOnManifold` to apply rotational diffusion to the
    particle quaternions :math:`\mathbf{q}_i` in simulations with active forces.
    The persistence length of an active particle's path is :math:`v_0 / D_r`.

Both clearly mention the relationship to the orientation quaternion q_i.