[lammps-users] Computing temperature gradient in a nanochannel filled with diatomic gas molecules

Dear Lammps Users,

I want to compute the temperature distribution in a nanochannel filled with diatomic gas molecules. My system is equilibrated at 300 K. When I use “compute temp” or “compute temp/com”, it indeed shows that the overall temperature in the gas domain is 300 K.
As you can see in the attached lammps input, I tried to derive the temperature distribution across the nanochannel using “compute chunk/atom”+ “fix ave/chunk” and also using “compute temp/chunk”, but using both approaches show that the temperature across the channel varies around 250 K.
The thing is previously I used the same approach for monoatomic gas and I haven’t faced such a problem in that case.
I highly appreciate it if someone could tell me how I can get the correct temperature distribution in a system with diatomic gas molecules.

#------------------------------------------input file--------------------------------------------

read_restart restart.Eq_after_thermos

group NiB region Ni_B

group NiT region Ni_T
#----------------------------------------------FORCE FIELD--------------------------------------------------------------

The discrepancy is caused by using fix shake, which removes degrees of freedom.
For the global temperature, that is easy to correct, since the number of removed degrees of freedom is known.
When using per-atom binning, however, this information is not readily available on a per-atom basis. How would you account for the loss of degrees of freedom, especially with complex molecules? so LAMMPS doesn’t even attempt to do that and disregards this. If you look at the documentation of compute temp/chunk, you see that you can set the number of degrees of freedom per atom manually. In your case that should be straightforward. You have two atoms and a bond between them, so you subtract half a DOF from each atom.