Problem in kinetic energy calculation

Dear LAMMPS users,

In order to equilibrate two flexible CH3 radical (with REBO potential), in a box with periodic boundary, at a given temperature 1100 K, I run LAMMPS with

fix 1 all nvt temp 1100.0 1100.0 0.1.

For output data, I use three groups of atoms (type 1 is C atoms and type 2 is H atoms):

group Catoms type 1
group Hatoms type 2
group CHatoms union Catoms Hatoms

So, the Catoms group contains NC=2 atoms and the Hatoms group contains NH=6 atoms.

and I compute separate temperatures for each of these groups:

compute TC Catoms temp
compute TH Hatoms temp
compute TCH CHatoms temp

In the same time I save in a dump file, the velocity components of all atoms.
When I check the total kinetic energy of the NC=2 C atoms, it corresponds exactly to 3/2*NC*k*TC, no problem.
BUT, when I check for H atoms, it does not correspond to 3/2*NH*k*TH, with NH=6 H atoms.
Moreover, and surprisingly (at least for me !), the H atoms total kinetic energy correspond exactly to 3/2*(NH-1)*k*TH, that is, as if it does not take into account the whole H atoms !
And it's the same for the CHatoms group, whose total kinetic energy is equal to 3/2*(NCH-1)*k*TCH.
Is there anyone who can help me understanding why this equivalence is ok for C atoms whereas it is not for H atoms and the whole atoms ?
Best regards.

Laurent.

Two questions

is this going to be a test system for a much larger calculation? If not you can easyly an ab initio model.

Have you checked the documentation on how translational degreesOf freedom ofThe system are handeled in your cases?

For 6 atoms the state of equilibrium is a bit fuzzy. Please consider a single harmonic oscillator. When is it in equilibrium?

Axel

Two questions

is this going to be a test system for a much larger calculation? If not you can easyly an ab initio model.

Yes.
This radical will then interact with a carbon substrate. It's an elementary process.
By doing a statistics of many such an elementary interaction, for different initial conditions of the CH3 radical (position, internal state) its sticking coefficient will be determined.
As the experiment is carried out at constant temperature (1100 K), both the radical and the substrate have to be equilibrated at 1100 K.
The equilibration step (at 1100 K) is to be followed by a production step to get an ensemble of initial states by sampling the NVT.

Have you checked the documentation on how translational degreesOf freedom ofThe system are handeled in your cases?

The assumption is made that the different degrees of freedom (trans, rot and vib) are in equilibrium.
From a previous posts,
http://sourceforge.net/mailarchive/forum.php?thread_name=CAENuAmFn2js7qW111WSs03u6ucnaKJvP0uH1BsmxDNt94MJ2fw%40mail.gmail.com&forum_name=lammps-users
I wall told that by using fix nvt I sould get good equipartition of energy into all the modes of energy of the CH3 radical.

From the lammps documentation:
The temperature is calculated by the formula KE = dim/2 N k T, where KE = total kinetic energy of the group of atoms (sum of 1/2 m v^2), dim = 2 or 3 = dimensionality of the simulation, N = number of atoms in the group, k = Boltzmann constant, and T = temperature.

But when comparing the total kinetic energy KE=sum of 1/2 m v^2 and comparing with 3/2NkT, for the three different groups (group of C atoms, group of H atoms and group of C and H atoms) I found discrepancy for H group and CH group. But not for C group.
But what really surprised me is that the equality
sum of 1/2 m v^2=3/2(N-1)kT
was fulfilled for H group and CH group instead of
sum of 1/2 m v^2=3/2N)kT !

For 6 atoms the state of equilibrium is a bit fuzzy. Please consider a single harmonic oscillator. When is it in equilibrium?

Axel
--
Dr. Axel Kohlmeyer [email protected] http://goo.gl/1wk0
International Centre for Theoretical Physics, Trieste. Italy.

Thanks for your reply.
Laurent.

If you want to see equipartioning with such a small system, you need a so-called massive thermostat, which is one where each degree of freedom is thermalized individually.
However in your scenario, this is in my opinion irrelevant. Sure the substrate should be well equilibrated, but in your experiment the radical is unlikely thermalized. What does it couple with? It will have some kinetic energy when generated and then likely have a high total energy in the vibrational degrees of freedom. Probably rotational, too. And keep in mind that lammps assumes that systems don't translate.

I sense you are trying to solve a problem that may not be as controllable as you assume.

Axel.

If you want to see equipartioning with such a small system, you need a so-called massive thermostat, which is one where each degree of freedom is thermalized individually.
However in your scenario, this is in my opinion irrelevant. Sure the substrate should be well equilibrated, but in your experiment the radical is unlikely thermalized. What does it couple with? It will have some kinetic energy when generated and then likely have a high total energy in the vibrational degrees of freedom. Probably rotational, too. And keep in mind that lammps assumes that systems don't translate.

Thanks a lot.
What do you mean by

"lammps assumes that systems don't translate"

?
Regards.
Laurent.

A periodic system is translationally invariant. Those are the three DOFs that are subtracted globally to get the temperature.

Axel