Effects and correction of angular momentum non-conservation in RNEMD for calculating thermal conductivity

Highlights

• Angular momentum flux in RNEMD may result in nonphysical wave motion heat transport.

• The wave leads to incorrect calculation of thermal conductivity at high heat flux.

• We improve RNEMD by choosing the hottest atoms in center-of-mass frames.

• The improved RNEMD eliminates wave motion and gives constant thermal conductivity.

Abstract

The effects of angular momentum non-conservation in reverse non-equilibrium molecular dynamics (RNEMD) for thermal conductivity calculation are investigated. The velocity swapping operation in RNEMD results in an abnormal angular momentum flux when the center-of-mass velocity of heat source is considerable. The angular momentum is found to motivate a nonphysical wave motion at high heat flux, which was regarded as a novel heat conduction mechanism in earlier researches. This additional heat transport channel may lead to a nonphysical increase in the computed thermal conductivity. Therefore, we propose an improved RNEMD, in which the hottest/coldest atoms are selected in center-of-mass frames. The results show that the improved RNEMD can greatly reduce the nonphysical correlation at any heat flux, eliminate the wave motion mode, and give constant thermal conductivity. Our research reveals a nonphysical effect from angular momentum non-conservation in RNEMD, and provides an improved RNEMD for heat transport in solid materials.

Introduction

Molecular dynamics (MD) simulation is a powerful way to study nanoscale heat transport. Based on Newton’s law of motion and interatomic potentials, it traces trajectories of all atoms, providing a direct approach to study nanoscale heat conduction mechanisms. Molecular dynamics has been successfully applied to study thermal conductivities [1], [2], [3], [4], phonon mean free path [5], phonon modes [6], [7] and heat transport across interfaces [8], [9], [10]. There are two typical methods to simulate thermal properties, i.e., Equilibrium Molecular Dynamics (EMD) and Non-Equilibrium Molecular Dynamics (NEMD). In EMD, the Green-Kubo method is commonly used to compute thermal conductivity from the heat current autocorrelation function. It is based on the linear response theorem [11], hence valid in the linear regime but invalid in the nonlinear regime [12]. In addition, long simulation time is needed to get converged heat current autocorrelation function [13], [14]. In NEMD, the thermal conductivity is computed directly from Fourier’s law

$$J_{\mu }=-\sum _{\nu }{k}_{\mu \nu }\frac{\partial T}{\partial x_{\nu }}.$$

This method requires a steady temperature gradient and heat current. One way is to establish the temperature gradient by keeping the heat source and heat sink at different temperatures by thermostats, like Berendsen [15], Nosé-Hoover [16], [17] and Langevin [18], [19] thermostats. Another typical way is Reverse Non-Equilibrium Molecular Dynamics (RNEMD) proposed by Müller-Plathe in 1997 [14], which introduces a heat flux between heat sink and heat source to motivate temperature gradient. It reverses the cause-and-effect relationship between heat flux and temperature gradient. The RNEMD method has been successfully used in various kinds of materials, such as silicon [20], grapheme [3], [21] and carbon nanotubes [1], [22], and was proved effective. Since the heat flux can be directly extracted from the algorithm, it greatly reduces the computation costs. Its convergence time is similar with other NEMD methods, much shorter than EMD.

The RNEMD method works as follows. For every several time steps, select $N$ hottest atoms in heat sink and $N$ coldest atoms in heat source, then swap their velocities. Since the velocity distributions in the slabs are wide enough, the hottest atoms in heat sink are expected to be hotter than the coldest atoms in heat source. This operation introduces a heat flux between the heat sink and the heat source, and finally reaches a stable state. By changing the frequency of swap operation and the number of atoms to be swapped, the heat flux applied to the system can be regulated. The heat flux in each swap operation is given by:

$$J=\frac{{\sum }_i^{N_{\mathit{swap}}}\frac{m}{2}\left(v_{\mathit{ih}}^2-v_{\mathit{ic}}^2\right)}{2t·L_y{L}_z},$$

where $N_{\mathit{swap}}$ is the number of atoms swapped, $m$ is the mass of atom, $v_{\mathit{ih}}$ is the velocity of $i\text{th}$ hottest atom in heat sink, $v_{\mathit{ic}}$ is the velocity of $i\text{th}$ coldest atom in heat source. $t$ is the time interval between each swap operation. $L_y$ and $L_z$ are the size of system along $y$ direction and $z$direction, respectively, assuming that the heat flux is along $x$ direction.

It is obvious that the total energy and momentum of the system are conserved in this algorithm, but the angular momentum is not conserved. Müller-Plathe guessed that this “should normally not be a problem” [14]. This is true since the angular momentum flux seems to be completely random in most cases. But no researches have investigated when the angular momentum flux will go nonrandom and its consequence. The influence of the angular momentum non-conversation in the simulations of heat conduction remains unclear.

Recently, using RNEMD, several researches reported a wave motion heat conduction mechanism at high heat flux in materials like carbon nanotubes [23], graphene [24], graphene nanoribbons [25] and hybrid graphene/silicene monolayers [26]. The wave motion transport mode is triggered at high heat flux and forms a wave-like temperature profile. It provides an additional channel of heat conduction and thus greatly increases thermal conductivity of these materials at high heat flux. In the view of phonons, a specific low frequency acoustic phonon mode is enhanced and goes ballistic, resulting in non-Fourier transport of heat. But unfortunately, this heat conduction mechanism was never demonstrated experimentally.

In this paper, we reveal that the wave motion heat conduction mode is just the consequence of the angular momentum non-conservation in RNEMD. The angular momentum introduced by RNEMD is found to regularly fluctuate at the frequency of the wave motion. By analyzing the Pearson correlation coefficient ($r$) between the angular momentum flux and the center-of-mass (COM) velocity of slabs, a strong linear correlation ($r>0.9$) between these two factors at high heat flux confirms that the angular momentum flux introduced by RNEMD is the source of the wave motion heat conduction. Further, we present a correction to the original RNEMD method to eliminate the regular angular momentum flux and keep the angular momentum flux random at any heat flux.