[lammps-users] compute heat/flux


I am not very familiar with statistical theories (and for that matter even with C++). So I have some basic questions and I’d really appreciate some clarification on these:

Considering the dimensions of the kappa equation: kappa= thermal conductivity, V = volume, kB=Boltzmann constant, T=temperature, t= time, J should be in a unit like W/m2 (or the corresponding unit system) which is the normal unit for heat flux. However, the equation J = Seivi + …… where ei is the per atom energy, vi is velocity, J has units of W-m or units of a flux-volume quantity.

  1. So, is the volume term in the kappa equation already absorbed in the J term as defined in the second equation?

For the input script given in the documentation, T= 70, V = 10213.3 (after the 8000 steps of fix npt), kB= 0.0019872067. Then, V/kBT2 = 1048.88. The integral of correlation term from the python script is 1.23e-10. Dt = 200 (fs). So the product of the 3 terms is 2.58e-5. On multiplying by the conversion factor (2917703220.0), it becomes 75283.9. This is ~3*105 times the value of 0.25 W/m-K given in the documentation.

  1. What am I missing in the above?

On the conversion factor itself:

With the following real units:

V: A3, kB : kcal/mol-K, T: K, t: fs, and J= kcal/mol-A2-fs (assuming flux quantity, rather than flux volume)

Unit of kappa : kcal/mol-K-fs-A = 4186 J/(6.023 x 1023 x 10-15 s x 10-10 m-K) = 69500 W/m-K. So the conversion factor from real to SI units should be ~69500.

  1. Could you please point out the error in this unit conversion and how the value listed in the documentation is obtained?

Thanks a lot to all and excuse me if I missed something(s) very simple (which is quite likely).


I would send your questions to Reese Jones (rjones at sandia.gov).


Dear Rohit,

I do not know about your calculations but here is mine. I have not calculated the integration of the heat flux vector but regarding the overall prefactor in front of the integration. Here it is.

From real to SI units, (kcal/mol-Angs-fs-K) to W/m-K, it is factor of 69467 (as you stated). Now you divide it with VolTemp^2kb (in kcal/mol-K). This comes to ~0.7. (Assuming volume is what you mentioned in your email). Then you have to multiply it with 40 (timestep of J output). This results in ~28. Now you have to multiply the overall intergration of J.J by this number. I am not fully sure if your integration value of ~1.e-10 is correct.

Reese: Any comments?


Just now, I looked at the script and it seems like what you are autocorrelating is J/vol not J. (That is why your integration is so low). In that case, if you multiply autocorrelation by (vol)^2. You get thermal conductivity about 0.6

Basically, looking at my last calculation.
28 (prefactor) * 1.23e-10 (integration) * (1.0213e+4)^2 (vol^2) = 0.59 W m-K. (This number will depend upon the quality of your integration).

I hope this will solve your concern.


Dear Vikas,

Thanks for both your responses. After your first response, I was
wondering why you've used the term V*T^2*kb instead of V/(T^2*kb). I
think your second email resolves other questions as well. I think you
were using a different approach to explain the calculation than that
given in the documentation.
However, I was able to get the no.s in the documentation after getting
cues from your email and following the documentation.

1. The problem was that I was initially interpreting the documentation
as 2917703220.0 being the conversion factor. After your response, it
seems that the term stands for the product of the three terms: the
Vol/(T^2*kb), the time step, and the conversion factor from
(kcal/mol-Angs-fs-K) to W/m-K. If this is correct, it would be great if
the documentation could be modified so that it is clearer to new users.
Assuming the above, I took another stab at the numbers:
Vol/(T^2*kb) ~ 1048.88, conversion factor ~ 69500, and time step = 40.0,
=> product of the 3 terms = 291588640 which is roughly the factor given
in the documentation.

Multiplying this with 1.23 x 10-10 gives thermal conductivity of roughly
0.359 W/m-K which is different from the value given in documentation,
but similar to your calculations (you did a calculation error in the
final step or probably a typo, missing the 3).

2. It seems that the J in the two equations in the documentation are
different quantities: J in the first equation is a flux quantity
(dimensions of W/m2), while the one in the 2nd equation is a flux-volume
(dimensions of W-m) quantity. If this is true, I think the two J's
should be denoted with different symbols, possibly J1 and J2 (or the
second equation should be J1*vol = ...)

Since the documentation explicitly asks to multiply the volume term
later to the integral of autocorrelation of flux, it seems to me that
the J given by the compute heat/flux should be a flux quantity
(dimensions of W/m2). However, looking at compute_heat_flux.cpp it seems
that the calculation is based on the 2nd equation of the documentation
and has not been normalized by the volume. That is, it is a flux-volume
quantity. Could you please clarify this apparent discrepancy?

I am not familiar with C++ and it is very likely that I am missing
something obvious. Please excuse me if that is the case.

Thanks again for your help.