Uncertainty quantification in non-equilibrium molecular dynamics simulations of thermal transport

Bulk thermal conductivity estimates based on predictions from non-equilibrium molecular dynamics (NEMD) using the so-called direct method are known to be severely under-predicted since finite simulation length-scales are unable to mimic bulk transport. Moreover, subjecting the system to a temperature gradient by means of thermostatting tends to impact phonon transport adversely. Additionally, NEMD predictions are tightly coupled with the choice of the inter-atomic potential and the underlying values associated with its parameters. In the case of silicon (Si), nominal estimates of the Stillinger-Weber (SW) potential parameters are largely based on a constrained regression approach aimed at agreement with experimental data while ensuring structural stability. However, this approach has its shortcomings and it may not be ideal to use the same set of parameters to study a wide variety of Si-based systems subjected to different thermodynamic conditions. In this study, NEMD simulations are performed on a Si bar to investigate the impact of bar-length, and the applied temperature gradient on the discrepancy between predictions and the available measurement for bulk thermal conductivity at 300 K by constructing statistical response surfaces at different temperatures. The approach helps quantify the discrepancy, observed to be largely dependent on the system-size, with minimal computational effort. A computationally efficient approach based on derivative-based sensitivity measures to construct a reduced-order polynomial chaos surrogate for NEMD predictions is also presented. The surrogate is used to perform parametric sensitivity analysis, forward propagation of the uncertainty, and calibration of the important SW potential parameters in a Bayesian setting. It is found that only two (out of seven) parameters contribute significantly to the uncertainty in bulk thermal conductivity estimates for Si.