Challenges in molecular dynamics simulations of heat exchange statistics

We study heat exchange in temperature-biased metal-molecule-metal molecular junctions by employing the molecular dynamics simulator LAMMPS. Generating the nonequilibrium steady state with Langevin thermostats at the boundaries of the junction, we show that the average heat current across a gold-alkanedithiol-gold nanojunction behaves physically, with the thermal conductance value matching the literature. In contrast, the full probability distribution function for heat exchange, as generated by the simulator, violates the fundamental fluctuation symmetry for entropy production. We trace this failure back to the implementation of the thermostats and the expression used to calculate the heat exchange. To rectify this issue and produce the correct statistics, we introduce single-atom thermostats as an alternative to conventional many-atom thermostats. Once averaging heat exchange over the hot and cold thermostats, this approach successfully generates the correct probability distribution function, which we use to study the behavior of both the average heat current and its noise. We further examine the thermodynamic uncertainty relation in the molecular junction and show that it holds, albeit demonstrating nontrivial trends. Our study points to the need to carefully implement nonequilibrium molecular dynamics solvers in atomistic simulation software tools for future investigations of noise phenomena in thermal transport.