Accelerating first-principles molecular-dynamics thermal conductivity calculations for complex systems

TL;DR

The paper tackles the computational bottlenecks of first-principles, Green-Kubo thermal-conductivity calculations in complex quasi-1D systems by developing a transferable MACE interatomic potential for InAs and implementing an efficient, unfolded heat-flux GK workflow. It systematically evaluates noise-reduction techniques, showing cepstral analysis accelerates convergence for low-conductivity ZB nanowires but underestimates high-conductivity WZ systems, where uncertainty propagation and covariance-aware methods prove essential. The authors introduce an uncertainty-based framework (KUTE) and an Euler-covariance approach to quantify and propagate HFACF uncertainties, enabling robust, quantitative error estimates across independent GK runs. The integrated workflow—combining MLIP-based MD, gauge-fixed heat flux, cepstral analysis (with model averaging), and covariance-aware uncertainty analysis—delivers a practical, accelerated, and robust path to accurate thermal-conductivity predictions across diverse nano-scale materials and structures.

Abstract

Atomistic simulations of heat transport in complex materials are costly and hard to converge. This has led to the development of several noise-reduction techniques applicable to equilibrium molecular-dynamics (MD) simulations. We analyze the performance of those strategies, taking InAs nanowires as our benchmark due to the diverse structures and complex phonon spectra of these quasi-1D systems. We demonstrate how, for low-thermal-conductivity systems, cepstral analysis can reduce computational demands while still delivering accurate results that do not require discarding arbitrary parts of the dataset. However, issues with this approach are revealed when treating high-thermal-conductivity systems, where the thermal conductivity is significantly underestimated. We discuss alternative methods to be used in that situation, relying on uncertainty propagation from independent simulations. We show that the contributions of the covariance matrix have to be included for a quantitative assessment of the error. The combination of these strategies with machine-learning interatomic potentials (MLIPs) provides an accelerated, robust workflow applicable to a diverse set of systems, as our examples using a highly transferable MACE potential illustrate.

Figures (10)

• Figure 1: Atomistic structures of the tested zincblende $[*]{111}$ and wurtzite $[*]{001}$ nanowires from different perspectives.
• Figure 2: Validation of the approaches for data generation and model training. (a) Correlation between mean force error and aggregated mean force uncertainty on the entire generated dataset, evaluated based on an ensemble of 5 MACE models that were trained on a fraction (4384 configurations) of said dataset. (b) Parity plot for the forces on all structures in the dataset, computed using the final model employed in the production MD simulations. (c) Phonon band structure for the InAs WZ phase as evaluated with the final MACE model and DFT. (d) Parity plot for a set of randomly displaced structures for the ZB $[*]{111}$ and WZ $[*]{001}$ nanowires.
• Figure 3: Comparison of different noise-reduction techniques for the ZB nanowire: using the full flux as in equation \ref{['main:e_hflux']} (blue line), removing the convective term leaving only the potential flux (red line), exploiting the gauge invariance of the flux as in equation \ref{['main:e_gauge']} (orange line), and employing a smoothening on the gauge fixed flux based on the lowest frequency in the vibrational density of states (green). Comparisons shown are of (a) the heat flux autocorrelation functions, (b) the thermal conductivity $\kappa$, (c) the power spectrum, where the relevant low time or frequency data are shown dominantly and the full ranges are shown as insets.
• Figure 4: Illustration of the GK evaluation procedure for the ZB nanowire. (a) Low-frequency power spectrum of the HFACF scaled to thermal conductivity units. The red curve indicates the effective smoothing of the spectrum when evaluating the thermal conductivity using the objectively chosen optimal number of cepstral coefficients $P$. (b) Second-order Akaike information criterion (AIC$_c$) as a function of the number of cepstral coefficients. The red dashed lines indicate the minimum. (c) Akaike weight distribution for the range of cepstral coefficients to be used for model averaging. (d) Thermal conductivity $\kappa$ as a function of the number of cepstral coefficients. The resulting values using the optimal number of cepstral coefficients and model averaging are indicated with the red and green dashed lines. (e) Time convergence curves for the thermal conductivity using the direct method and the cepstral analysis. The curves are given for different correlation lengths. (f) Cutoff dependence of the cepstral analysis using the full simulation time.
• Figure 5: Comparison of different noise-reduction techniques for the WZ nanowire: using the full flux as in equation \ref{['main:e_hflux']}, removing the convective term leaving only the potential flux, exploiting the gauge invariance of the flux as in equation \ref{['main:e_gauge']}, and for additionally applying a smoothening filter with a window width corresponding to the lowest peak in the vibrational density of states. Comparisons shown are for (a) the heat flux autocorrelation functions, and (b) the thermal conductivity $\kappa$.