[lammps-users] Reg: Calculating translational and rotational energy per bin of nitrogen molecule in channel flow

Dear All,
I am simulating flow of nitrogen gas molecules bounded by two walls. I was wondering is there a command or a way to calculate translational and rotational energies of molecules separately per bin (The channel is divided into many 1D bins). I only found the command bond/local which calculates energies for the entire system, but I want to calculate per bin.

It would be of great help if you could help figure out a method to calculate the per bin energies.

k kishore kumar

You would have to do some C++ programming and implement suitable per-chunk computes that compute rotational and translational kinetic energy per chunk.
Even then the collection into bins would be a bit convoluted, since you first need to define chunks per molecule and then compute those two terms based on those chunks.
Then you would need to use compute spread/chunk/atom to redistribute those values back to the atoms and then use a second compute chunk/atom to do the 1d binning and collect the per-chunk energies for rotation and translation. Please note that with this approach you would get double the energies because you are copying the per-molecule energies to both constituent atoms and thus need to scale the results by 0.5.

Since you have rather simple molecules, it may be simpler to do this analysis in post-processing from a trajectory that has atom positions, atom and molecule IDs and velocities. Then you can easily compute the per molecule total kinetic energy and the kinetic energy of the center of mass and by subtracting that from the total kinetic energy get the rotational kinetic energy component.


Thanks a lot for the reply, Axel.

The second suggestion seems to be easy to implement. I wrote a code to calculate both translational and rotational energies per bin. Surprisingly, the sum (translational energy of molecule + rotational energy as calculated below) doesn’t add up to the total energy calculated from considering thermal velocities of each atom.
The way I am calculating the rotational energy is as follows (logic taken from bond/local)
inertia[i,j] += delxdelx + delydely + delzdelz
omegasq[i,j] = (delcvx
delcvx + delcvydelcvy + delcvzdelcvz)/(delxdelx + delydely + delz*delz)

ike_rot[i,j] = 0.5*inertia[i,j]*omegasq[i,j]

delx = x[i,k] - x_mol[i,j] (length of bond)
delcvx = cu[i,k] - cvx_mol[i,j] (difference between thermal velocity of atom and center of mass thermal velocity of molecule)

Could you please tell me if there is something wrong with this logic?

Thanks once again in advance.

k kishore kumar

The code you refer to computes 3 components: COM translation, vibration, and rotation.

I have no time to debug your code fragments.

I don’t want to get this types of massage please remove me from cc

You have subscribed yourself to the mailing list and you have to unsubscribe yourself.
The link to the mailing list management interface is given in the footer of every email.

My email archive also shows that you have already been given the shortcut link to the page where you can directly unsubscribe (or change delivery settings). You just have to do it. Nobody will do it for you. Sourceforge’s privacy policy does not permit mailing list admins to see the list of subscribers or change their settings there.