Temperature from vx, vy, vz profiles (ave/chunk)

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

Nevermind, I figured it out. For future reference (for newcomers):

vx vy vz in ave/chunk are averaging over chunk atoms as well as over
time. The average of values distributed around zero either cancel out
completely, or result in very small values (pretty much what should
happen in an equilibrated system).

Therefore dumping vx vy vz in ave/chunk will be useless if we wish to
calculate temperature from it. At the minimum, we need to square them
before any averaging happens. Therefore we should dump vxsq (vx
squared) and vysq, vzsq instead. We can do it like so:

compute vel all property/atom vx vy vz
variable vxsq atom c_vel[1]*c_vel[1]
variable vysq atom c_vel[2]*c_vel[2]
variable vzsq atom c_vel[3]*c_vel[3]

compute 6 all chunk/atom bin/1d z lower 1.0 units box
fix 7 all ave/chunk 1 10000 20000 6 v_vxsq v_vysq v_vzsq file out.dat

And now we can use the squared velocity components outside lammps to
calculate temperature, and this should give the correct result.

Thank you and best regards,
F