I have a question that probably is very silly, and might not only limit to lammps but to all MD calculations in general.

I am trying to calculate the kinetic energy and heat capacity for an ensemble of gas molecules and I find the the way MD calculations produces the thermodynamic quantities not quite adequate. In lammps, the kinetic energy of the ensemble is basically calculated by summing the 1/2mv^2 of all atoms, this will result in kinetic energy equals to 3/2NRT, N been the number of atoms in the ensemble, and result in constant volume heat capacity Cv equals to 3/2R normalized to per atoms. However, this strategy would fail if my ensemble is made by non monoatomic molecules.

For example, for diatomic molecules, the ideal Cv at ambient temperature would be 5/2R (translational energy + rotational energy) and 7/2R at high temperature (vibrational energy kicks in). Lammps MD calculation would simply give 3R if normalized to per molecules.
For mon-linear triatomic molecules, the ideal Cv at ambient temperature would be 3R (translational energy + rotational energy) and 6R at hight temperature (vibrational energy kicks in). Lammps MD calculation would simply give 4.5R if normalized to per molecules.

The point here is is there a way to perform lammps MD simulations based on dynamics of "molecules" instead of "atoms", and extract the information of rotational and vibrational energies from such simulations?

The point here is is there a way to perform lammps MD simulations based on dynamics of "molecules" instead of "atoms", and extract > the information of rotational and vibrational energies from such simulations?

The dynamics of molecules with flexible bonds are no different
than the dynamics of atoms. The temp of the collection
of atoms is equal to the temp of the collection of molecules.
You could post process a dump file and compute vibrational
and rotational energy if you want. They should be equal to
the temperature for an equilibrated system.

I shouldnt comment because im not an expert in this field , and probably i will make the ridiculous, but anyways here it goes.

Many thermodynamics properties can be obtained from the “density of states” (DOS) as example in the Einstein model the DOS is a delta function and At high temperatue limit we obtain Cv = 3NKb. Thus From the density of states we can obtain the vibrational helmholtz energy or the vibration entropy and thus derive the specific heat (Cv).

We can calculate the density of states from Molecular dynamics simulations (MD) and there are different methods in the literathure one of such methods is using the velocity autocorrelation function.

Thank you very much for the insight Steve! I do believe the dynamics of individual atoms shall be same as the flexible bond molecules.

But I am getting this concern when I see the behavior of the constant heat capacity Cv changes in the fashion directly proportional to the atom numbers of the system. For instance, at ambient temperature when vibrational energy contribution to the heat capacity is little, an ideal gas system composed by N non-linear molecules that are made by 10 atoms each molecule should have ideal Cv equals to 3R, contributed from 3 translational and 3 rotational degrees of freedom for each molecule. But lammps MD simulation would give Cv equals to 15R because there are total 10N atoms and each atom contribution 3/2R to Cv.

I didn't run a long simulation for a molecule that big but I did tried some calculations up to tetra-atomic molecules in lammps with the reaxff force filed at a short period of time, they all suffer the same issue I am seeing. I also explicitly did the MD simulation of O2 molecule , using 0.1 fs time step running for 1 million steps, and I get Cv equals to 3R (normalized to number of molecules) instead of 5/2R for such a diatomic linear molecule, and I also don't Cv been excited to 7/2R or any where close at high temperature (3000K), it stays at 3R all the time. Is there a post processing approach I need to take to get this sorted out?

Thanks for the suggestion on DOS. I calculated Cv in two different methods. The first is taking the gradient of the internal energy against temperature. Another method is using the Fluctuation-dissipation theorem, calculating the Cv from the energy fluctuation. The two methods give the consistent results that both suffer the issue I discussed previously.

Hi Qing,
I never did this type of calculation but I think there might be a flaw
in your reasoning. In real life due to the quantum nature of the
system certain modes will not be active at low temperatures thus no
energy can be allocated into these nodes at all. You need energy jumps
to get there as the density of states is discontinuous. Nothing like
his happens in classical MD where you can always "excite" a vibration
mode no matter how small is the energy injected in the system. To make
the comparison fair I think you have to freeze your bonds in the MD
simulation and introduce the corresponding number of constraints in
the definition of the heat capacity. This is what my intuition tells
me here. Not sure though.
Carlos