Hello all,
I think I have found a potential error in fix_rigid, in the langevin thermostat which is built into that fix.
As detailed in the documentation, the langevin option adds two terms to the force and torque on each rigid body: a random force/torque and a drag term proportional to the velocity of the body. This has the effect of thermostating the rigid bodies, and can be used e.g. to perform langevin/brownian dynamics. For the rotational motion the term added to the torque about each of x,y,z axis (in FixRigid::post_force() in fix_rigid.cpp) is:
langextra[i][3] = inertia[i][0]*gamma1*omega[i][0] + sqrt(inertia[i][0])*gamma2* (random->uniform()-0.5);
The problem arises because, as I understand it, lammps stores the principal moments of inertia of each rigid body in the array inertia. There are three values for each body - the moments of inertia about each of the principal axis of the body, i.e. the axes which constitute the frame of reference of the body (so the array inertia does not change with time). When the orientation of the body is integrated, then e.g. to calculate the angular momentum from the angular velocity, the angular velocity must first be found in the frame of reference of the body. That is to say, the array omega and the torques are in the frame of reference of the system, whereas inertia is in the frame of reference of the body.
So should, in FixRigid::post_force(), omega first be converted to the body frame of reference, then the langextra force calculated in the body frame, before being converted back into the system frame of reference? I've tried to do some tests of this, which are detailed in the attached pdf; I also edited FixRigid::post_force() to conform to the above, and compared results with those of the original version.
Also, one further suggestion about the thermostating of the motion of the rigid bodies. A recent update (24 Apr 2013) to the fix langevin for the case of extended particles allows the user to specify that the damping time for thermostating of translational and rotational velocities are different (by specifying the factor by which they are different). This can be useful when using this fix to do langevin dynamics. It would also be useful to have this option in the fix rigid langevin option.
Regards,
Chris Brackley
tests.pdf (74.8 KB)
Hmm - the lingo LAMMPS uses is body-frame
and space-frame (Cartesian axes of the simulation box).
Please take a look at these two methods,
which should be doing the same thing, one for
single ellipsoidal particles, one for rigid bodies
FixLangevin::angmom_thermostat()
FixRigid::post_force()
In FixLangevin, only the space-frame angmom
and quat are stored, so a space-frame omega
is derived from them via the call to mq_to_omega()
(also see the code for that in math_extra.cpp)
In FixRigid, omega is already the space-frame
omega, so this transformation isn’t needed.
I thougt is was correct to apply the new Langevin
torque in the space frame, isn’t it? I.e. the
rotational temperature of the particle or body
is computed from the space-frame omega
and torques should always be in the space frame.
But you are saying this is incorrect, i.e. the
thermostat should be applied in the body-frame?
Steve
Hi Steve,
Thanks for your reply.
I think that you are right, it is fine to apply the Langevin torque in the space-frame. However, the problem is in that the drag term contains the moment of inertia - and I think this is in the body-frame.
In FixRigid::post_force() for example there is a term:
langextra[i][3] = inertia[i][0]*gamma1*omega[i][0] + sqrt(inertia[i][0])*gamma2*(random->uniform()-0.5);
and I think this could be incorrect since inertia[i][0] is the moment of inertia about one of the principal axes of the ith body, i.e. in the body-frame. Whereas omega[i][0] is the rotational velocity about one of the axes in the space-frame. I may be wrong, but think for this term to make sense, everything must be in the same frame. So I guess the correct way would be to first calculate the moment of inertia about each of the axes in the space frame (and this will change as the body rotates), and then calculate langextra[i][3] etc. (I think in my original post I described doing things the other way round, i.e. transforming omega to body-frame, calculate langextra in body-frame, then transform back to space-frame; I think that amounts to the same thing, but is probably a bit less efficient, and not so clear.)
I'm not so familiar with FixLangevin::angmom_thermostat(), but it could be that it is already taken care of there, as I see that the inertia array is used in calculating omega, which already involves a change of frame.
Chris
I may be wrong, but think for this term to make sense, everything must be in the same frame.
I think this is right - I wasn’t thinking about the inertia terms.
FixLangevin is doing the same thing, b/c this eq:
inertia[0] = EINERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
is simply the moment of inertia for an ellipsoid around one of its principal axes.
So the 2 classes are doing the same thing, but both could be wrong (when
omega and angmom are not aligned).
Let me think about this a bit more …
Steve