I have a question about the equation to calculate thermal conductivity using the outputted information from “fix thermal/conductivity” command.
According to the manual, “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 appropriate units”.
However, the paper by Muller-Plathe mentioned the heat flux should be half of the energy transferred by these swaps divided by the cross section because the heat could flow back to two directions from the hottest slab equally. I am wondering whether I should divide the value from “fix thermal/conductivity” by 2 to get the thermal flux or the LAMMPS has already accounted for the factor internally.
Please help me verify whether the equation I derived is correct. Thank you.
Units= metal
Simulation box =8X8X8
Lattice parameter 5.47
So cross section = (8*5.47)^2
CET = cumulative energy transfer outputted by “fix thermal/conductivity”
Nevery is the parameter in “fix thermal/conductivity”
Ntotal=# of swaps conducted by “fix thermal/conductivity”
2= factor accounts for equal probability of heat conduction in two direction
Looking at the source code, I don’t see that factor of 2 you mention, so you’ll probably need to include it in your post-processing analysis in order to follow the prescription in the paper you cite. You can confirm this by carefully examining the source code and the literature.
I think that you haven’t got that thermal conductivity formula quite right. Looks to me like you should divide the total energy exchanged by the total time elapsed, regardless of the swapping interval or rate. Adapting from your formula, I’d favor something more like this:
The document for the fix states that…The scalar is the cummulative kinetic energy transferred between the bottom and middle of the simulation box (in the edim direction) is stored as a scalar quantity by this fix
I think we do not need to divide by 2 because the energy transferred between one end and the middle (hot) layer is only stored in the fix, unlike the total energy transferred throughout the entire domain when the factor of 2 might come into play. Please correct me if I am wrong.