# Rotational kinetic energy

I am a new Lammps user
I am running a simulation of ellipsoids (hybrid bond ellipsoid), with an orientation-dependent potential that I created.
In the example here, there is no translational move of the particles, but there is an orientational move, and this should produce a rotational kinetic energy. However, in the thermo logfile it remains 0.

The way I do it is :
compute 1 all erotate/asphere
thermo_style custom step ke pe ebond etotal temp

and the thermo logfile is like this
Step KinEng PotEng E_bond TotEng Temp
0 0 5 5 5 0
10 0 4.9989633 4.9989633 4.9989633 0
20 0 4.9834332 4.9834332 4.9834332 0
30 0 4.9165807 4.9165807 4.9165807 0
40 0 4.7401431 4.7401431 4.7401431 0
50 0 4.3845898 4.3845898 4.3845898 0
60 0 3.7920856 3.7920856 3.7920856 0
70 0 2.9535454 2.9535454 2.9535454 0

The orientation-dependent potential energy seems OK. Also, the angular momentum, torque and quaternion are nonzero in the dumps:

ITEM: ATOMS id type x y z ix iy iz xu yu zu angmomx angmomy angmomz tqx tqy tqz quatw quati quatj quatk
1 1 0 0 -1 0 0 0 0 0 -1 0 0 1.97543 0 0 4.38633 0 0.780806 0.624773 0
2 1 0 0 0 0 0 0 0 0 0 0 0 -1.97543 0 0 -4.38633 0 0.624773 0.780806 0

What is incorrect in the way I compute the kinetic energy ?
The ellipsoids have a shape 111 (spheres, but I want to store the orientation). I do not need to define a tensor of inertia ?

Thanks !

Sam

The default compute that thermo output uses for KE
is compute temp, which is just translational eng.
You can create your own compute with compute temp/asphere
or compute erotate/asphere and print those values with
thermo output. And you can use thermo_modify to replace
the default output temp with one of those.

Steve

Thank you for your help. I was able to compute the rotational part.

Now I would like to study the Brownian dynamics of my system. This is usually done by fix langevin + fix nve. But I want to include the rotational degrees of freedom.
I have two questions/problems:

1. I see that in the new version of Langevin, there is an option (angmom) to include the rotations in the thermostat.
However, in the doc, I don't see that there is a corresponding rotational damping. Does this damping exist ?

2. I built a new orientation-dependent bond, and in the previous version of Lammps (where Langevin operates only on translations) I could compile without a problem. However, with the new version, the definitions of aspherical particles (ellipsoids in my case) seem to be different. The mathematical functions on quaternions (stored in math_extra.h) have different names, which is not a big problem. But when I try to get the quaternion of the atom, I get the following message:
bond_ejtehadi.cpp:73: error: 'class LAMMPS_NS::Atom' has no member named 'quat'
The corresponding line in the bond file is :
double **quat = atom->quat;
and it seems that there is no "quat" parameter in this new version. What should I change ?

Thanks

Sam

Thank you for your help. I was able to compute the rotational part.

Now I would like to study the Brownian dynamics of my system. This is
usually done by fix langevin + fix nve. But I want to include the
rotational degrees of freedom.
I have two questions/problems:

1. I see that in the new version of Langevin, there is an option
(angmom) to include the rotations in the thermostat.
However, in the doc, I don't see that there is a corresponding
rotational damping. Does this damping exist ?

2. I built a new orientation-dependent bond, and in the previous version
of Lammps (where Langevin operates only on translations) I could compile
without a problem. However, with the new version, the definitions of
aspherical particles (ellipsoids in my case) seem to be different. The
mathematical functions on quaternions (stored in math_extra.h) have
different names, which is not a big problem. But when I try to get the
quaternion of the atom, I get the following message:
bond_ejtehadi.cpp:73: error: 'class LAMMPS_NS::Atom' has no member named
'quat'
The corresponding line in the bond file is :
double **quat = atom->quat;
and it seems that there is no "quat" parameter in this new version. What
should I change ?

have a look at, e.g. pair_gayberne.cpp
the information you are looking for is accessed
there via:

AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid;

and

iquat = bonus[ellipsoid[i]].quat;

HTH,
axel.

Comment below.

Steve

Thank you for your help. I was able to compute the rotational part.

Now I would like to study the Brownian dynamics of my system. This is
usually done by fix langevin + fix nve. But I want to include the
rotational degrees of freedom.
I have two questions/problems

1. I see that in the new version of Langevin, there is an option
(angmom) to include the rotations in the thermostat.
However, in the doc, I don't see that there is a corresponding
rotational damping. Does this damping exist ?

See the method angmom_thermostat in fix_langevin.cpp.
It works exactly the same as the point particle thermostat,
adding the 2 gamma terms to the torque on the particle.
So that includes the damping, I believe.

Thank you for your help. I was able to compute the rotational part.

Now I would like to study the Brownian dynamics of my system. This is
usually done by fix langevin + fix nve. But I want to include the
rotational degrees of freedom.
I have two questions/problems:

1. I see that in the new version of Langevin, there is an option
(angmom) to include the rotations in the thermostat.
However, in the doc, I don't see that there is a corresponding
rotational damping. Does this damping exist ?

2. I built a new orientation-dependent bond, and in the previous version
of Lammps (where Langevin operates only on translations) I could compile
without a problem. However, with the new version, the definitions of
aspherical particles (ellipsoids in my case) seem to be different. The
mathematical functions on quaternions (stored in math_extra.h) have
different names, which is not a big problem. But when I try to get the
quaternion of the atom, I get the following message:
bond_ejtehadi.cpp:73: error: 'class LAMMPS_NS::Atom' has no member named
'quat'
The corresponding line in the bond file is :
double **quat = atom->quat;
and it seems that there is no "quat" parameter in this new version. What
should I change ?

have a look at, e.g. pair_gayberne.cpp
the information you are looking for is accessed
there via:

AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid;

and

iquat = bonus[ellipsoid[i]].quat;

Thanks, I managed to compile
What about the rotational damping for Langevin ?
The translational part may lead to a reasonable dynamics in the case where the rotational and translational dof are strongly coupled, but for a general use it would be good to have a rotational part too, isn't it ?

From the doc of the fix langevin command, I understand that this

orientational part is not implemented. Can you confirm this or tell me how to use it ?
Sam

Assuming you saw my answer from a few minutes ago,
I said I think damping is implemeted for rotation (both the
Ff and Fr terms in the Langevin eq on the doc page).

But maybe I don't understand your Q. Why do you think it is not?

Steve

Well, the doc says
The keyword angmom and omega keywords enable thermostatting of rotational degrees of freedom in addition to the usual translational degrees of freedom.
I interpreted the term “thermostatting” as referring to the Fr random force only, but I admit it may be my mistake.
I will check the code of the fix to be sure
Thanks

Sam

I have a problem to dump the orientation of the particles with the new version of the ellipsoids
Before, I used dump custom quatw quatx ...
Now it seems that the quaternion is no longer an option for dump custom
How can I get the orientation of the particles ?
Sam

no - it's both terms. Most people (including me) think
of Langevin as a thermostat. That's all it means.

But check the code to be sure it's what you want.

Steve

Compute property/atom covers every pre-atom property
that LAMMPS has. Dump custom is just the more common
ones, b/c you can always include a per-atom computes properties in dump custom.

Steve