Heat flux profile calculations

Snippets for heat flux profile calculations

This snippet illustrates the implementation of heat flux calculations using LAMMPS, centroid.

Please refer to the article: Polarization of ethanol under thermal fields, Guansen Zhao, Yaxuan Xiao, Angad Deshmukh, Aidan Chapman and Fernando Bresme J. Chem. Phys. 163, 024501 (2025)

--- Heat flux calculations ---

Note

Note that the approach below uses the so-called IK1 route. This is an approximation. However, it has been shown to produce heat flux profiles in agreement with the fluxes obtained from the continuity equation. See e.g. Water Polarization under Thermal Gradients. For systems involving large density gradients, e.g., in fluids under nanoconfinement conditions one should consider the Method of planes.

The local heat flux within a small control volume $V$ can be calculated using,

$$\mathbf{J}_q = \frac{1}{V} \left( \sum_{i \in V} e_i \mathbf{v}_i - \sum_{i \in V} \boldsymbol{\sigma}_{i} \cdot \mathbf{v}_i \right),$$

where,

• $i\in V$ are the atoms within the control volume, $V$.
• $e_i$ is the energy of atom $i$.
• $\boldsymbol{\sigma}_i$ is the per atom stress tensor.
• $\mathbf{v}_i$ is the velocity of atom $i$.

Stress decomposition (real-space + k-space contributions)

The compute centroid/stress/atom does not support the calculation of kspace contributions. Support can either be patched in or by using a workaround by creating an additional stress/atom compute that only calculates the kspace contributions:

compute stress_virial all centroid/stress/atomNULL pair bond angle dihedral improper fix
compute stress_kspace all stress/atom NULL kspace

Per-atom energies

These are simply calculated with builtin computes in lammps.

compute myKE all ke/atom
compute myPE all pe/atom

Combining the components

Conversion of stress: atm·A^3 -> kcal/mol (units real)

variable atmAAcubed2kcal equal 1.458e-5

Heat current (IK1): e.g. x component: $J_x = \sum_i [ (e_i) v_{x,i} - (σ_{xx} v_{x,i} + σ_{xy,i} v_{y,i} + σ_{xz,i} v_{z,i}) ]$

e_i = KE_i + PE_i [kcal/mol] σ_* from compute stress/atom virial, units atm·Å^3; convert to kcal/mol

variable Jx_a atom (c_myKE+c_myPE)*vx-(c_stress_virial[1]*vx+c_stress_virial[4]*vy+c_stress_virial[5]*vz)*v_atmAAcubed2kcal
variable Jy_a atom (c_myKE+c_myPE)*vy-(c_stress_virial[4]*vx+c_stress_virial[2]*vy+c_stress_virial[6]*vz)*v_atmAAcubed2kcal
variable Jz_a atom (c_myKE+c_myPE)*vz-(c_stress_virial[5]*vx+c_stress_virial[6]*vy+c_stress_virial[3]*vz)*v_atmAAcubed2kcal

variable Jx_k atom -(c_stress_kspace[1]*vx+c_stress_kspace[4]*vy+c_stress_kspace[5]*vz)*v_atmAAcubed2kcal
variable Jy_k atom -(c_stress_kspace[4]*vx+c_stress_kspace[2]*vy+c_stress_kspace[6]*vz)*v_atmAAcubed2kcal
variable Jz_k atom -(c_stress_kspace[5]*vx+c_stress_kspace[6]*vy+c_stress_kspace[3]*vz)*v_atmAAcubed2kcal

Chunking along z (heat flux direction)

variable bw equal $(lz/300) # Here assuming 300 bins
compute zchunk_all all chunk/atom bin/1d z 0 ${bw} units box

fix chunkave_flux all ave/chunk 1 $(v_nsteps) $(v_nsteps) zchunk_all v_Jx_a v_Jy_a v_Jz_a v_Jx_k v_Jy_k v_Jz_k file ${outputname}.flux_profiles ave running

Post-simulation Analysis

From the calculation above we get $J_x$ [(Kcal/mol)(Å/fs) "units real"]

conversion-factor = Convert (kcal/mol)*Å/fs per atom -> W (J m/s)

nbins -- number of z bins (match your zchunk_all)

rho_x -- number density in the bin [particles/m3]

Sum kinetic + potential (per-chunk averages) -> multiply by Ncount to get total per chunk and divide by volume_bin to get the flux in [W/m^2]

$J_{x} = (J_{x,a} + J_{x,k}) * \text{conversion-factor} * \rho_x$