# Tutorial: Thermal transport from NEMD and HNEMD

## Introduction

• In this tutorial, we use the NEMD and HNEMD methods to study heat transport in graphene at 300 K and zero pressure, from the ballistic to the diffusive transport regime. The spectral decomposition method as described in Ref. [1] is used here.
• All the input and output files can be found here.

## Preparing the Inputs

### The inputs for the NEMD simulation

#### The xyz.in file

• We consider a graphene sheet of size of about 25 nm x 43 nm (40400 atoms). The transport is in the $y$ direction. We divide the length in the $y$ direction into 9 groups. Group 0 (a small one) will be fixed. Group 1 (about 9 nm long) will act as a heat source and group 8 (about 9 nm long) will act as a sink. The remaining middle part is evenly divided into a few groups (group 2 to group 7).

#### The run.in file

• The run.in input file is given below:
potential    potentials/tersoff/Graphene_Lindsay_2010_modified.txt 0
velocity     300

# equilibration
ensemble     nvt_ber 300 300 0.01
fix          0
time_step    1
dump_thermo  1000
run          100000

# production
ensemble     heat_lan 300 100 10 1 8
fix          0
compute      0 10 100 temperature
compute_shc  4 2 250 1
run          1000000

• The second line of command with the velocity keyword tells that the velocities will be initialized with a temperature of 300 K.
• There are two runs. The first run serves as the equilibration stage.
• Here, the NVT ensemble (the Berendsen thermostat) is used. The target temperature is 300 K and the coupling constant is 0.01 (dimensionless) for the thermostat.
• The time step for integration is 1 fs.
• Atoms in group 0 will be fixed.
• The thermodynamic quantities will be output every 1000 steps.
• There are $10^5$ steps (100 ps) for this run.
• The second run is for production.
• Here, two local Langevin thermostats are applied to generate a nonequilibrium heat current. The time parameter in the Langevin thermostats is 0.1 ps (100 time steps). The target temperature of the heat source is $300+10=310$ K and the target temperature of the heat sink is $300-10=290$ K. Atoms in group 1 are in the heat source region and atoms in group 8 are in the heat sink region.
• Atoms in group 0 will be fixed.
• The compute keyword is used to compute the group temperatures and energy transfer rate. Here, grouping method 0 is used, and the relevant data are sampled every 10 time steps and averaged for every 100 data points before written out.
• The line with the compute_shc keyword is used to compute the spectral heat current (SHC). The SHC will be calculated for group 4 in the default grouping method 0. The relevant data will be sampled every 2 steps and the maximum correlation time is $250 \times 2 \times 1 {\rm fs} = 500 {\rm fs}$. The transport directions is 1 (y direction).
• There are $10^6$ steps (1 ns) in the production run. This is just an example. To get more accurate results, I suggest you use $10^7$ steps (10 ns).

### Inputs for the HNEMD simulation

#### The xyz.in file

• The same as in the NEMD simulation, but no atoms will be fixed in the run.in file.

#### The run.in file

The run.in file reads:

potential    potentials/tersoff/Graphene_Lindsay_2010_modified.txt 0
velocity     300

# equilibration
ensemble     nvt_nhc 300 300 100
time_step    1
dump_thermo  1000
run          1000000

# production
ensemble      nvt_nhc 300 300 100
compute_hnemd 1000 0 0.00001 0
compute_shc   4 2 250 1
run           10000000

• The first two lines are the same as in the NEMD simulation.
• There are two runs. The first run serves as the equilibration stage.
• Here, the NVT ensemble (the Nose-Hoover chain thermostat) is used. The target temperature is 300 K and the thermostat coupling constant is 0.1 ps (100 time steps).
• The time step for integration is 1 fs.
• The thermodynamic quantities will be output every 1000 steps.
• There are $10^6$ steps (1 ns) for this run.
• The second run is for production.
• Here, the global temperature is controlled by the Nose-Hoover chain thermostat with the same parameters as in the equilibration stage.
• The compute_hnemd keyword is used to add a driving force and compute the thermal conductivity using the HNEMD method [1]. Here, the conductivity data will be averaged for each 1000 steps before written out, and the driving force parameter is $10^{-5}$/A and is in the $y$ direction.
• The line with the compute_shc keyword is used to compute the SHC. The relevant data will be calculated for group 4 in the default grouping method 0, sampled every 2 steps, with a maximum correlation time of $250 \times 2 \times 1 {\rm fs} = 500 {\rm fs}$, and the transport direction is 1 (y direction).
• There are $10^7$ steps (10 ns) in the production run. This is just an example. To get more accurate results, I suggest you use $2 \times 10^7$ steps (20 ns) and do a few independent runs and then average the relevant data.

## Results and discussion

### Computation time

• The NEMD simulation takes about 15 min (Tesla K40) or a few min (Tesla P100).
• The HNEMD simulation takes about 3 hours (Tesla K40) or about 1 hour (Tesla P100).

### NEMD Results

• TODO: Figures to be updated.
(a) Temperature profile in the NEMD simulation. (b) Energies accumulated in the thermostats.
• The figure above shows the temperature profile and the energies accumulated in the thermostats. The energy of the thermostat coupling to the heat source region is decreasing, because energy is transferred from the thermostat to the atoms in the source region. The energy of the thermostat coupling to the heat sink region is increasing, because energy is transferred from the atoms in the sink region to the thermostat. The absolute values of the slopes of the lines in (b) should be the same; otherwise it means energy is not conserved. The absolute slope is the energy transfer rate $Q=dE/dt$.
• The thermal conductance in this system can be calculated as

$$G = \frac{Q/S}{\Delta T}$$ where $S$ is the cross-sectional area and $\Delta T$ is the temperature difference, which is 19 K here (slightly smaller than the target value of 20 K). The calculated thermal conductance is about 10 GW m-2 K-1. This is the classical value. For quantum correction, see Ref. [2].

(a) Force-velocity correlation function. (b) Spectral thermal conductance.
• The figure above shows the force-velocity correlation function $K(t)$ and the spectral thermal conductance $G(\omega)$. For theoretical background, see the spectral decomposition method.

### HNEMD results

• Add more results here.

(a) Force-velocity correlation function $K(t)$. (b) Spectral thermal conductivity $\kappa(\omega)$. (c) Spectral phonon mean free path $\lambda(\omega)$. (d) Length-dependent thermal conductivity $\kappa(L)$.
• The above figure shows the results from the HNEMD simulation [1].
• (a) The force-velocity correlation function $K(t)$.
• (b) The spectral thermal conductivity $\kappa(\omega)$.
• (c) The spectral phonon mean free path calculated as $\lambda(\omega)=\kappa(\omega)/G(\omega)$, where $G(\omega)$ is the quasi-ballistic spectral thermal conductance from the above NEMD simulation.
• (d) The length-dependent thermal conductivity calculated as

$$\kappa(L) = \int \frac{d\omega}{2\pi} \frac{\kappa(\omega)}{1 + \lambda(\omega)/L}.$$

## References

• [1] Zheyong Fan, Haikuan Dong, Ari Harju, and Tapio Ala-Nissila, Homogeneous nonequilibrium molecular dynamics method for heat transport and spectral decomposition with many-body potentials, Phys. Rev. B 99, 064308 (2019). https://doi.org/10.1103/PhysRevB.99.064308
• [2] Zhen Li, Shiyun Xiong, Charles Sievers, Yue Hu, Zheyong Fan, Ning Wei, Hua Bao, Shunda Chen, Davide Donadio, and Tapio Ala-Nissila, Influence of Thermostatting on Nonequilibrium Molecular Dynamics Simulations of Heat Conduction in Solids, J. Chem. Phys. 151, 234105 (2019). https://doi.org/10.1063/1.5132543