Setting drag coefficients for rigid bodies with Langevin dynamics

[crisbutu]

Hello,

I have a question related to handling drag for rigid bodies when using Langevin dynamics.

1. Is the drag coefficient of a rigid body determined solely by the gamma parameter of its center bead, ignoring constituent beads?

2. I'm uncertain about the interpretation of the gamma_r parameter (rotational drag coefficient). The documentation for version 2.9.7, which I'm using, specifies its units as mass/time (the same for v3). With the default setting being gamma_r = (1,1,1), does this imply that the rotational diffusivity of the rigid body is equivalent to that of a sphere, regardless of its actual shape? Or does the internal handling consider the geometry of the rigid body to produce an expected drag tensor with units of mass*length^2/time?

Specifically, I want to implement rigid rods and have set the moment of inertia for the central bead as that of a cylinder. I want to ensure that the rotational diffusivity about the axial direction differs from the perpendicular directions, as expected for a rod. Could anyone advise on how to correctly configure this?

Any help would be greatly appreciated. Thank you!

[joaander]

The three components of gamma_r are the rotational drag coefficients about the three axes of rotation (x, y, z). gamma_r has units of 1/time because each component appears in the equation of motion as:
$$I \frac{d\vec{L}}{dt} = \vec{\tau}_C - \gamma_r \cdot \vec{L} + \vec{\tau}_R$$
(this equation is in the documenation). Note that \cdot here is an element-wise multiplication - not a dot product.

The angular momentum $L$ has units [mass][length]^2[time]^(-1) and torque $\tau$ has units [mass][length]^2[time]^(-2), so gamma_r must have units of [time]^(-1).

However, looking at the code - the documentation for this may be incorrect:

hoomd-blue/hoomd/md/TwoStepLangevin.cc

Lines 316 to 355 in 103d1b4

	// original Gaussian random torque
	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.0, 0.0);
	Scalar rand_x = hoomd::NormalDistribution<Scalar>(sigma_r.x)(rng);
	Scalar rand_y = hoomd::NormalDistribution<Scalar>(sigma_r.y)(rng);
	Scalar rand_z = hoomd::NormalDistribution<Scalar>(sigma_r.z)(rng);
	// check for degenerate moment of inertia
	bool x_zero, y_zero, z_zero;
	x_zero = (I.x == 0);
	y_zero = (I.y == 0);
	z_zero = (I.z == 0);
	bf_torque.x = rand_x - gamma_r.x * (s.x / I.x);
	bf_torque.y = rand_y - gamma_r.y * (s.y / I.y);
	bf_torque.z = rand_z - gamma_r.z * (s.z / I.z);
	// ignore torque component along an axis for which the moment of inertia zero
	if (x_zero)
	bf_torque.x = 0;
	if (y_zero)
	bf_torque.y = 0;
	if (z_zero)
	bf_torque.z = 0;
	// change to lab frame and update the net torque
	bf_torque = rotate(q, bf_torque);
	h_net_torque.data[j].x += bf_torque.x;
	h_net_torque.data[j].y += bf_torque.y;
	h_net_torque.data[j].z += bf_torque.z;
	if (D < 3)
	h_net_torque.data[j].x = 0;
	if (D < 3)
	h_net_torque.data[j].y = 0;

I don't have time right now to fix it. Look at rand_x (which is a function gamma_r.x) and gamma_r.x * (s.x / I.x) and figure out what units it should be. Feel free to file a pull request that makes the correction to the docs.

[zmcdargh]

I think this equation is a mistake, and the L here is supposed to be angular velocity, not angular momentum (the net torque on a body is just dL/dT, not I*dL/dt). Making this change resolves the issue with the dimensions of gamma_r, the discrepancy between the code and the docs, and a few inconsistencies in the docs. Either I or @crisbutu will do a PR to fix this. Thanks for your help!

[joaander]

Thanks! I look forward to the PR.