I'm trying to understand temperature calculation from velocities (for

my understanding of how lammps does it).

I'm running a vapor-liquid-vapor LJ system (2.5 cutoff). It's in

equilibrium but thermostat is applied only to bulk of the liquid (not

the interface, and not to vapor). It's in 3D. N is about 125,000, and

volume is about (55, 55, 223) in sigma. (constant volume sim).

I decided to calculate temperature myself (outside of lammps) by

dumping velocity profiles (w.r.t z-axis) using ave/chunk. So I used

the following commands:

compute 6 all chunk/atom bin/1d z lower 1.0 units box

fix 7 all ave/chunk 1 10000 20000 6 vx vy vz file out.dat

Now since 3kT/2 = m * average(vx^2 + vy^2 + vz^2) / 2

For LJ, k=1, and I set m=1. Thermostat is at T=0.8. And after

cancellation, rearrangement, we get:

T = average(vx^2 + vy^2 + vz^2) / 3

However when I calculate right-hand-side (RHS), it is about 1500-2000

times smaller, i.e., I'm getting values of 0.0005 for RHS, when T is

0.8 (reduced LJ units).

I plotted the profiles, w.r.t z, and the middle (liquid part) has

practically zero velocity. But even after ignoring that, T from vapor

velocities still comes out 1000 times smaller.

Could you help me with the calculation? Is the big difference due to

"bias" since I've read we need to remove bias from velocities to

accurately calculate T. Or am I missing something serious (like some

conversion factor, or taking density or number of particles into

account)?

Thanks in advance,

F