Hi everybody,
I am student and not that experienced with LAMMPS, so sorry if this question is silly or I missed something in the documentation.
I am doing a simulation of Hexanitrostilbene with the ReaxFF force field, based on the files that are given in the LAMMPS examples of reaxff/HNS. My first goal now is to measure the vibrational energy of the molecules, but I’m not sure how to do that. I found some paper that say I can calculate it by calculating the sum of the kinetic energy of the atoms in the molecule and then subtracting the kinetic energy of the center of mass and the rotational energy. But I’m not completely sure how to calculate the moment of inertia for the rotational energy and so I’m wondering if there is a simpler, “direct” way of calculating the vibrational energy from simulation outputs?
Greetings and thanks in advance!
Florian
To do this processing on-the-fly with LAMMPS is going to be a bit tricky.
First of all, with ReaxFF you have no explicit molecules since the bond topology can change when reactions happen (which is the main point where ReaxFF is different from other empirical force fields). So you need to use fix reaxff/species to partition the system into “molecules” (according to the computed bond order) and you can use the molecule IDs from the per-atom property of the fix in combination with compute chunk/atoms to have the system partitioned accordingly.
Once you have the chunk information, you can run various per-chunk computes to pre-compute several useful per-chunk properties and write them to a dump file, as well as use compute ke/atom to get the per-atom kinetic energy and then you can post-process the whole lot.
However, it may be easier to just have one custom dump file that contains the (unwrapped!) positions, velocities, and the molecule partitioning from fix reaxff/species and then you would do the processing by yourself in a custom script/program.
Hello Axel,
thank you for your reply!
Currently, I am doing it mostly like you said, I am using the reaxff/bonds, evaluate the molecules in my own script and then calculating all the properties I need from the dumped positions and velocities.
So it sounds like I can stick to this way of calculation.
Greetings
Florian