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] = (delcvxdelcvx + 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.