Exact heat flux formula and its spectral decomposition in molecular dynamics for arbitrary many-body potentials

TL;DR

This work derives an exact heat-flux framework for molecular dynamics with arbitrary many-body potentials and provides a spectral decomposition (SDHC) to dissect mode contributions. Grounded in Hardy and Torii formalisms with energy partitioning via weights, the method computes both control-volume and control-surface heat currents and their spectral content, applicable to diverse materials. Validation against Green-Kubo and Non-Equilibrium MD shows that the centroid-based heat flux yields accurate thermal conductivities, while conventional LAMMPS formulations significantly underestimate them for several layered and covalent systems. The SDHC results reveal pronounced size-dependent evolution of phonon contributions, notably enhanced ZA-mode transport in graphene at larger lengths and dominant acoustic contributions in MoS$_2$, underscoring the importance of exact formalisms for phonon-mediated heat transfer in complex potentials.

Abstract

In this study we have derived an exact framework for the calculation of the heat flux and its spectral decomposition in Molecular Dynamics (MD) for arbitrary many-body potentials. This work addresses several lacks and limitations of previous approaches and allows for the accurate computational study of thermal properties in a wide variety of many-body systems with MD. We have tested our modifications with Green-Kubo (GK) and Non-Equilibrium MD (NEMD) simulations for various 2D and 3D material systems using the Tersoff and Stillinger-Weber potentials as examples. The spectral decomposition of the heat current was also calculated for monolayer graphene (1LG) and MoS2 , for different system lengths. Our results show that the heat current calculated by our method is consistently in agreement with the thermostat current in NEMD, while previous implementations can estimate quite poorly the thermal conductivity both under GK and NEMD simulations, and both for 2D and 3D materials. The decomposition of the heat current also sheds light on the contribution of different phonon modes to thermal conductivity and its dependence on length. Our methodology is implemented in the widely used LAMMPS code specifically for the Tersoff and SW potentials, and it is readily applicable to the vast majority of many-body MD potentials.

Figures (5)

• Figure 1: Visualization of the HF calculation workflow for a 3-body interaction with arbitrary weights $p_i$. (a) Average HF in a Control Volume: the 3-body potential energy $U$ is distributed to $i$, $j$, $k$ according to their weights $p_i$ (left) and the atomic stress tensors are calculated (middle) according to the 'old' group formula in eq. (\ref{['eq:LAMMPS_Stress_Tensor']}) (gray) and the centroid formula from eqs. (\ref{['eq:Centroid_Stress_Tensor']}) and (\ref{['eq:Centroid_Position']}) (black). The HF is calculated in a control volume for both cases from eq. (\ref{['eq:HeatFlux_virial']}) (right) (b) HF across a Control Surface: $U$ is distributed to regions $1$ and $2$ according to the weights $p_i$ of the atoms belonging to them (left) and the forces $\boldsymbol{F}_i^{1 \leftarrow 2}$ are calculated (middle) accordingly from eq. (\ref{['eq:F12']}) (black). No previous such implementation exists except for pair forces (gray). $Q^{\text{virial}}_{1 \rightarrow 2}$ is then calculated from eq. (\ref{['eq:HeatCurrent_Torii']}) (right), from which other quantities can be derived, such as $Q(\omega)$. Although bulk a-SiO$_2$ was chosen for visualization, the methodology is applicable to any arbitrary system.
• Figure 2: Comparative results with the Green-Kubo method. The calculated thermal conductivity $\kappa$ of various 2D and 3D systems using GK and $\boldsymbol{J}$ calculated from the centroid formula of eq. (\ref{['eq:Centroid_Stress_Tensor']}) (in pink) and the group formula of eq. (\ref{['eq:LAMMPS_Stress_Tensor']}) (in blue). $\kappa_{group}$ is shown as a percentage of $\kappa_{centr}$. For the centroid formula, the absolute values of $\kappa$ (in $W/mK$) are also displayed. In all cases, the uncertainties were of the order of 2-5%
• Figure 3: Comparative results with the NEMD method. The heat currents from the mean heat flux values calculated with the centroid formula of eq. (\ref{['eq:Centroid_Stress_Tensor']}) ($Q_{centr}$, in pink) and the group formula of eq. (\ref{['eq:LAMMPS_Stress_Tensor']}) ($Q_{group}$, in blue), shown as a percentage of the current $Q_{thermo}$ calculated from the energy exchange rate of the thermostats. For the values of $Q$ from the HF, the uncertainties were of the order of 2-5%, while for $Q_{thermo}$ it was below 1%.
• Figure 4: The spectral decomposition $Q(\omega)$ for 1LG for different system lengths. (a) $Q(\omega)$ for 1LG of length $L$=50 nm (in black). The in-plane component $Q^{\text{in}}$ (in red) and the out-of-plane component $Q^{\text{out}}$ (in blue) are also shown. (b) $Q^{\text{in}}$ (in red) and $Q^{\text{out}}$ (in blue) components are shown for different system lengths: $L$=50 nm (solid lines) and $L$=500 nm (dashed lines). (c) The phonon DOS of 1LG (in black), along with the in-plane (in red) and out-of-plane (in blue) contribution, normalized to unit total area.
• Figure 5: The spectral decomposition of the heat current $Q(\omega)$ for MoS$_2$ with different system lengths (a) $Q(\omega)$ for $L$=50 nm (red) and $L$=500 nm (black). (c) The total phonon DOS of MoS$_2$, normalized to unit total area.