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
Thanks in advance,