Problems with compute com (centre of mass velocity)

Dear LAMMPS users and developers,

I am having problems with the command “compute com” in a simulation of a lj system. I just want to calculate the translational and non-translational contributions to the total kinetic energy of all the atoms in the box.
When I do this modifying the LAMMPS code (in FixHeat::end_of_step(), for test purposes), I get a sensible value as compared to the total kinetic energy.

However, when I calculate the com velocity using “compute all com” and subsequently square and multiply it by the total mass over two, I get a value that is several orders of magnitudes off (117355.16).
I saw in the documentation of temp “compute com” that the vector is given in distance units. I am working with a lj system, that’s why several conversion factors are unity. Division by the timestep (0.001) would also not account for that.

My question is now: Why are the velocities calculated with group->vcm(0,group->mass(0), all_vcm) and compute all com not identical (group-id 0 is “all”)?

My apologies for posting such a simple question, but I got stuck.

Best wishes

Peter

Here is the (quick & dirty) piece of code I added to “fix_heat” (with which I am experimenting around) to create the additional output for validation.

double all_masstotal = group->mass(0);
double all_ke = group->ke(0)*force->ftm2v;
double all_vcm[3];
group->vcm(0,all_masstotal,all_vcm);

double all_vcm2 = (all_vcm[0]*all_vcm[0] + all_vcm[1]*all_vcm[1] + all_vcm[2]*all_vcm[2]);

double all_ke_trans = 0.5all_masstotalall_vcm2 * force->ftm2v;
double all_ke_ntrans = all_ke - all_ke_trans;

if (update->ntimestep % 200 == 0) {
printf("(FIXHEAT: all_vcm^2 = %lf, \t all_vcm = (%lf, %lf, %lf))\n",all_vcm2, all_vcm[0], all_vcm[1], all_vcm[2]);

}

Here is the console output (please ignore the first part in which the crystal is melted):

pw359@…4757…:~/fixheat/lj/test/300_dt0.001$ ~/apps/lammps-25Nov13_test/src/lmp_openmpi -in lj.lmp
LAMMPS (25 Nov 2013)
Lattice spacing in x,y,z = 1.6796 1.6796 1.6796
Created orthogonal box = (0 0 0) to (10.0776 10.0776 20.1552)
1 by 1 by 1 MPI processor grid
Created 1728 atoms
Setting up run …
Memory usage per processor = 2.76758 Mbytes
Step Temp KinEng PotEng TotEng Press
0 3 4.4973958 -6.9999116 -2.5025157 -4.2703406
200 1.6845059 2.5252967 -4.9835801 -2.4582834 5.2630082
400 1.6907258 2.534621 -4.9418267 -2.4072057 5.4643114
600 1.6763134 2.513015 -4.9063603 -2.3933453 5.626113
800 1.7144289 2.5701551 -4.9305688 -2.3604136 5.5360477
1000 1.7305419 2.5943106 -4.940948 -2.3466374 5.4973514
Loop time of 4.85131 on 1 procs for 1000 steps with 1728 atoms

Pair time () = 4.22079 (87.0031) Neigh time () = 0.464731 (9.57951)
Comm time () = 0.0492308 (1.01479) Outpt time () = 0.000480175 (0.00989784)
Other time (%) = 0.116078 (2.39272)

Nlocal: 1728 ave 1728 max 1728 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 6791 ave 6791 max 6791 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 242270 ave 242270 max 242270 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 242270
Ave neighs/atom = 140.203
Neighbor list builds = 50
Dangerous builds = 0
WARNING: New thermo_style command, previous thermo_modify settings will be lost (…/output.cpp:665)
Setting up run …
Memory usage per processor = 2.77502 Mbytes
Step KinEng ke_all1 ke_all_t vcm2_all vcm_all[ vcm_all[ vcm_all[
1000 2.5943106 2.5943106 117355.16 135.82773 4.5944417 4.6220746 9.6620529
(FIXHEAT: all_vcm^2 = 0.001971, all_vcm = (-0.043345, 0.001141, 0.009516))
(FIXHEAT: all_vcm^2 = 0.001971, all_vcm = (-0.043345, 0.001141, 0.009516))
1200 2.5751388 2.5751388 117320 135.78704 4.5857728 4.6223029 9.663956
(FIXHEAT: all_vcm^2 = 0.001971, all_vcm = (-0.043345, 0.001141, 0.009516))
(FIXHEAT: all_vcm^2 = 0.001971, all_vcm = (-0.043345, 0.001141, 0.009516))
1400 2.5802637 2.5802637 117284.98 135.74651 4.5771039 4.6225312 9.6658591
(FIXHEAT: all_vcm^2 = 0.001971, all_vcm = (-0.043345, 0.001141, 0.009516))
(FIXHEAT: all_vcm^2 = 0.001971, all_vcm = (-0.043345, 0.001141, 0.009516))
1600 2.5686618 2.5686618 117250.1 135.70613 4.568435 4.6227595 9.6677623

Relevant parts of the Input file (see full file attached lj.lmp)

fix NVE all nve
compute ke all ke/atom

1) Calculate total kinetic energy

compute ke_all1 all ke

2) Alternative way to calculate total kinetic energy (for validation)

compute ke_all2 all reduce sum c_ke

3) Alternative way to calculate the total kinetic energy (validation for 4))

variable ke3 atom mass*(vxvx+vyvy+vz*vz)/2.
compute ke_all3 all reduce sum v_ke3

4) Calcualte translational kinetic energy of entire box (ke_all_t)

compute vcm_all all com
compute m all property/atom mass
compute m_all all reduce sum c_m
variable vcm2_all equal c_vcm_all[1]*c_vcm_all[1]+c_vcm_all[2]*c_vcm_all[2]+c_vcm_all[3]c_vcm_all[3]
variable ke_all_t equal c_m_all
v_vcm2_all/2.
variable ke_all_nt equal c_ke_all1-v_ke_all_t

thermo 200
thermo_style custom step ke c_ke_all1 v_ke_all_t v_vcm2_all c_vcm_all[1] c_vcm_all[2] c_vcm_all[3]
thermo_modify flush yes
run ${Nprod}

lj.lmp (2.96 KB)

Dear LAMMPS users and developers,

I am having problems with the command "compute com" in a simulation of a lj
system. I just want to calculate the translational and non-translational
contributions to the total kinetic energy of all the atoms in the box.
When I do this modifying the LAMMPS code (in FixHeat::end_of_step(), for
test purposes), I get a sensible value as compared to the total kinetic
energy.

However, when I calculate the com velocity using "compute all com" and
subsequently square and multiply it by the total mass over two, I get a
value that is several orders of magnitudes off (117355.16).
I saw in the documentation of temp "compute com" that the vector is given in
distance units. I am working with a lj system, that's why several conversion
factors are unity. Division by the timestep (0.001) would also not account
for that.

My question is now: Why are the velocities calculated with
group->vcm(0,group->mass(0), all_vcm) and compute all com not identical
(group-id 0 is "all")?

My apologies for posting such a simple question, but I got stuck.

i am confused by your question. you are comparing the center of mass
(compute com) to the center of mass *velocity* group->vcm().

why _should_ those be identical. those are different properties...

you need to compare to group->xcm()

axel.

Hi Axel,

of course, you are right! For some reason I thought that compute com would compute the com-velocity. Even after re-reading the documentation twice. I am sorry for that!

Thanks,
Peter

Dear Axel,

I changed my input script now to calculate the com-velocity and the translational kinetic energy and non-trans. kinetic energy.
The com-velocity is now the same I get with group->vcm(), which is good.

However, the way I calculate translational and non-translational kinetic energy is still wrong. Something is wrong with the way I calculate the total mass, I think.
Am I using invalid combinations of computes and variables that are evaluated at different times?

Best wishes,

Peter

lj.lmp (1.81 KB)

Dear Axel,

I changed my input script now to calculate the com-velocity and the
translational kinetic energy and non-trans. kinetic energy.
The com-velocity is now the same I get with group->vcm(), which is good.

However, the way I calculate translational and non-translational kinetic
energy is still wrong. Something is wrong with the way I calculate the total
mass, I think.

i don't think so. you are probably confused by reduced units
defaulting to thermo_modify norm yes

axel.

Thanks again Axel! That did the trick!

Best wishes,
Peter