# question on " fix thermal/conductivity" command

Dear Lammps gurus:

Lammps provides a convenient tool for computing thermal conductivity via the “fix thermal/conductivity” command. The idea is to create a temperature gradient across simulation box by effectively imposing a heat flux across box (achieved via kinetic energy swapping).

The manual says " the total kinetic energy transferred by these swaps is computed by the fix and can be output. Dividing this quantity by time and the cross-sectional area of the simulation box yields a heat flux. The ratio of heat flux to the slope of the temperature profile is the thermal conductivity of the fluid, in appopriate units"

These sentences suggest that k = E_swap/(A dT/dx), where Eswap is the energy transfer rate via swapping. However, since the total heat flux is equally split between the two halfs of the simulation box, should we use 1/2E_swap for the heat flux instead? Or, the E_swap reported by Lammps already takes into account the 1/2 factor?

Thank you very much.

With best regards.

Nick

no - you need to put in the 1/2 (for a periodic system). Please
refer to the Muller-Plathe paper cited.

Steve

Hi Steve and all,

I found the recent discussion about 1/2 factor in fix_thermal/conductivity. I also used this command to calculate TC of graphene and other materials. Today I double-checked the subroutine of " fix_thermal_conductivity.cpp", if we search " eswap", we will find the exchange energy change actually has been divided by the factor of 1/2.

After swapping, the new velocity for site 1 and 2 are (2vm-v1) and (2vm-v2), respectively, so the kinetic energy change for cold/hot region is

dE=1/2m1[(2vm-v1)^2-v1^2]=1/2m1*(2vm-2v1)(2vm)=2m1(vm-v1)*vm.

But in the subroutine, the eswap is calculated in a form of m1*(vm-v1)*vm.

So I argue the 1/2 factor has been divided in the swap energy calcualtion. Therefore the heat flux may be calculated as eswap/(t. A), where t is the time and A is the cross-section area.

Is it correct?

Best!

Dongshan

There are several factors of 2 in this calculation, and the formulas
are more complicated with the extension to allow swaps between

I thought the factor of 2 you were asking about was in the denominator
of eq 4 in the original 1997 Muller-Plathe paper. That factor
of 2 is due to a periodic system (the energy diffuses both directions).

The numerator should simply be the KE that is swapped, i.e, the
KE added to one of the slabs by a velocity exchange.

If 2 particles with velocity 8 and 10 (in one dir) and mass=1 are swapped, then
the delta KE of one particle (slab) is 1/2 m 10^2 - 1/2 m 8^2 = 18.
The eswap formula in the code is delta KE = m * (9 * (9-8)), where vcm = 9.
So that is delta KE = 9, but each of the 2 atoms contributes,
so that is delta KE = 18.

So it looks right to me. I.e. the code is tallying swapped KE. To compute
eq 4 in the paper, you need to add the final factor of 2 in the denominator,
assuming you a have a periodic system.

I'm CCing Craig (who generalized the fix to treat different masses), to
see if he agrees.

Steve

Hi Prof. Craig,

Could you please say something to this problem? If the 1/2 factor should be added in the calculation of the thermal conductivity?

Dongshan