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?

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

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.

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

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.