Heat transport¶

EMD method¶

A popular approach for computing the lattice thermal conductivity is to use equilibrium molecular dynamics (EMD) simulations and the Green-Kubo (GK) formula. In this method, the running thermal conductivity (RTC) along the \(x\)-direction (similar expressions apply to other directions) can be expressed as an integral of the heat current autocorrelation (HAC) function:

\[\kappa_{xx}(t) = \frac{1}{k_BT^2V} \int_0^{t} dt' \text{HAC}_{xx}(t').\]

Here, \(k_{\rm B}\) is Boltzmann’s constant, \(V\) is the volume of the simulated system, \(T\) is the absolute temperature, and \(t\) is the correlation time. The HAC is

\[\text{HAC}_{xx}(t)=\langle J_{x}(0)J_{x}(t)\rangle,\]

where \(J_{x}(0)\) and \(J_{x}(t)\) are the total heat current of the system at two time points separated by an interval of \(t\). The symbol \(\langle \rangle\) means that the quantity inside will be averaged over different time origins.

Related keyword in the run.in file: compute_hac
Related output file: hac.out
Related tutorial: Thermal transport from EMD

We only used the potential part of the heat current. If you are studying fluids, you need to output the heat currents (potential and kinetic part) using the compute keyword and calculated the HAC by yourself.

We have decomposed the potential part of the heat current into in-plane and out-of-plane components [Fan2017]. If you do not need this decomposition, you can simply sum up some components in the hac.out file.

NEMD method¶

Non-equilibrium molecular dynamics (NEMD) can be used to study thermal transport. In this method, two local thermostats at different temperatures are used to generate a non-equilibrium steady state with a constant heat flux.

If the temperature difference between the two thermostats is \(\Delta T\) and the heat flux is \(Q/S\), the thermal conductance \(G\) between the two thermostats can be calculated as

\[G = \frac{Q/S}{\Delta T}.\]

Here, \(Q\) is the energy transfer rate between the thermostat and the thermostated region and \(S\) is the cross-sectional area perpendicular to the transport direction.

We can also calculate an effective thermal conductivity (also called the apparent thermal conductivity) \(\kappa(L)\) for the finite system:

\[\kappa(L) = GL = \frac{Q/S}{\Delta T/L}.\]

where \(L\) is the length between the heat source and the heat sink. This is to say that the temperature gradient should be calculated as \(\Delta T/L\), rather than that extracted from the linear part of the temperature profile away from the local thermostats. This is an important conclusion in [Li2019].

To generate the non-equilibrium steady state, one can use a pair of local thermostats. Based on [Li2019], the Langevin thermostatting method is recommended. Therefore, the ensemble keyword with the first parameter of heat_lan should be used to generate the heat current.

The compute keyword should be used to compute the temperature profile and the heat transfer rate \(Q\).
Related output file: compute.out
Related tutorial: Thermal transport from NEMD and HNEMD

HNEMD method¶

The homogeneous non-equilibrium molecular dynamics (HNEMD) method for heat transport by Evans has been generalized to general many-body potentials [Fan2019]. This method is physically equivalent to the EMD method but can be computationally faster.

In this method, an external force of the form [Fan2019]

\[\boldsymbol{F}_{i}^{\rm ext} = E_i \boldsymbol{F}_{\rm e} + \sum_{j \neq i} \left(\frac{\partial U_j}{\partial \boldsymbol{r}_{ji}} \otimes \boldsymbol{r}_{ij}\right) \cdot \boldsymbol{F}_{\rm e}\]

is added to each atom \(i\), driving the system out of equilibrium. According to [Gabourie2021], it can also be written as

\[\boldsymbol{F}_{i}^{\rm ext} = E_i \boldsymbol{F}_{\rm e} + \boldsymbol{F}_{\rm e} \cdot \mathbf{W}_i\]

Here, \(E_i\) is the total energy of particle \(i\). \(U_i\) is the potential energy of particle \(i\). \(\mathbf{W}_i\) is the per-atom virial. \(\boldsymbol{r}_{ij}\equiv\boldsymbol{r}_{j}-\boldsymbol{r}_{i}\), and \(\boldsymbol{r}_i\) is the position of particle \(i\).

The parameter \(\boldsymbol{F}_{\rm e}\) is of the dimension of inverse length and should be small enough to keep the system within the linear response regime. The driving force will induce a non-equilibrium heat current \(\langle \boldsymbol{J} \rangle_{\rm ne}\) linearly related to \(\boldsymbol{F}_{\rm e}\):

\[\frac{\langle J^{\mu}(t)\rangle_{\rm ne}}{TV} = \sum_{\nu} \kappa^{\mu\nu} F^{\nu}_{\rm e},\]

where \(\kappa^{\mu\nu}\) is the thermal conductivity tensor, \(T\) is the system temperature, and \(V\) is the system volume.

A global thermostat should be applied to control the temperature of the system. For this, we recommend using the Nose-Hoover chain thermostat. So one should use the ensemble keyword with the first parameter of nvt_nhc.

The compute_hnemd keyword should be used to add the driving force and calculate the thermal conductivity.
The computed results are saved to the kappa.out file.
Related tutorial: Thermal transport from NEMD and HNEMD

Spectral heat current¶

In the framework of the NEMD and HNEMD methods, one can also calculate spectrally decomposed thermal conductivity (or conductance). In this method, one first calculates the following virial-velocity correlation function [Gabourie2021]:

\[\boldsymbol{K}(t) = \sum_{i} \left\langle \mathbf{W}_i(0) \cdot \boldsymbol{v}_i (t) \right\rangle,\]

which reduces to the non-equilibrium heat current when \(t=0\).

Then one can define the following Fourier transform pairs [Fan2017]:

\[\tilde{\boldsymbol{K}}(\omega) = \int_{-\infty}^{\infty} dt e^{i\omega t} K(t)\]

where

\[\boldsymbol{K}(t) = \int_{-\infty}^{\infty} \frac{d\omega}{2\pi} e^{-i\omega t} \tilde{\boldsymbol{K}}(\omega)\]

By setting \(t=0\) in the equation above, we can get the following spectral decomposition of the non-equilibrium heat current:

\[\boldsymbol{J} = \int_{0}^{\infty} \frac{d\omega}{2\pi} \left[2\tilde{\boldsymbol{K}}(\omega)\right].\]

From the spectral decomposition of the non-equilibrium heat current, one can deduce the spectrally decomposed thermal conductance in the NEMD method:

\[G(\omega) = \frac{2\tilde{\boldsymbol{K}}(\omega)}{V\Delta T}\]

with

\[G = \int_{0}^{\infty} \frac{d\omega}{2\pi} G(\omega).\]

where \(\Delta T\) is the temperature difference between the two thermostats and \(V\) is the volume of the considered system or subsystem.

One can also calculate the spectrally decomposed thermal conductivity in the HNEMD method:

\[\kappa(\omega) = \frac{2\tilde{\boldsymbol{K}}(\omega)}{VTF_{\rm e}}\]

with

\[\kappa = \int_{0}^{\infty} \frac{d\omega}{2\pi} \kappa(\omega).\]

where \(F_{\rm e}\) is the magnitude of the driving force parameter in the HNEMD method.

This calculation is invoked by the compute_shc keyword and the results are saved to the shc.out file.

Related tutorial: Thermal transport from NEMD and HNEMD

Modal analysis methods¶

A system with \(N\) atoms will have \(3N\) vibrational modes. Using lattice dynamics, the vibrational modes (or eigenmodes) of the system can be found. The heat flux can be decomposed into contributions from each vibrational mode and the thermal conductivity can be written in terms of those contributions [Lv2016]. To calculate the modal heat current in GPUMD, the velocities must first be decomposed into their modal contributions:

\[\boldsymbol{v}_i (t) = \frac{1}{\sqrt{m_i}} \sum_{n} \boldsymbol{e}_{i,n} \cdot \boldsymbol{\dot{X}}_n(t)\]

Here, \(\boldsymbol{\dot{X}}_n\) is the normal mode coordinates of the velocity of mode \(n\) \(m_i\) is the mass of atom \(i\) \(\boldsymbol{e}_{i,n}\) is the eigenvector that gives the magnitude and direction of mode \(n\) for atom \(i\) \(\boldsymbol{v}_i\) is the velocity of atom \(i\)

The heat current can be rewritten in terms of the modal velocity to be:

\[\boldsymbol{J}^{\text{pot}} = \sum_{i} \mathbf{W}_i \cdot \left[ \frac{1}{\sqrt{m_i}} \sum_{n} \boldsymbol{e}_{i,n} \cdot \boldsymbol{\dot{X}}_n(t) \right] = \sum_{n} \left(\sum_{i} \frac{1}{\sqrt{m_i}} \mathbf{W}_i \cdot \boldsymbol{e}_{i,n} \right) \cdot \boldsymbol{\dot{X}}_n(t)\]

This means that the modal heat current can be written as:

\[\boldsymbol{J}^{\text{pot}}_n = \left(\sum_{i} \frac{1}{\sqrt{m_i}} \mathbf{W}_i \cdot \boldsymbol{e}_{i,n} \right) \cdot \boldsymbol{\dot{X}}_n(t)\]

This modal heat current can be used to extend the capabilities of the EMD and HNEMD methods. The extended methods are called Green-Kubo modal analysis (GKMA) [Lv2016] and homogeneous non-equilibrium modal analysis (HNEMA) [Gabourie2021].

Green-Kubo modal analysis¶

The Green-Kubo Modal Analysis (GKMA) calculates the modal contributions to thermal conductivity by using [Lv2016] [Gabourie2021]:

\[\kappa_{xx,n}(t) = \frac{1}{k_BT^2V} \int_0^{t} dt' \langle J_{x,n}(t')J_{x}(0)\rangle.\]

Here, \(k_{\rm B}\) is Boltzmann’s constant, \(V\) is the volume of the simulated system, \(T\) is the absolute temperature, and \(t\) is the correlation time. \(J_{x}(0)\) is the total heat current and and:math:J_{x,n}(t’) is the mode-specific heat current of the system at two time points separated by an interval of \(t'\). The symbol \(\langle \rangle\) means that the quantity inside will be averaged over different time origins.

Related input file: eigenvector.in
Related keyword in the run.in file: compute_gkma
Related output file: heatmode.out

For the GKMA method, we only used the potential part of the heat current.

HNEMDEC method¶

A system with \(M\) components has \(M\) independent fluxes: the heat flux and any \(M-1\) momentum fluxes of the \(M\) component. The central idea of Evans-Cummings algorithm is designing such a driving force that produce a dissipative flux that is equivalent to heat flux or momentum flux. By measuring the heat current and momentum current, we obtain onsager coefficents that can be used to derive the thermal conductivity.

In the case of heat flux \(\boldsymbol{J}_{q}\) as dissipative flux, for the \(i\) atom belonging to \(\alpha\) component:

\[\boldsymbol{F}_{i}^{\alpha,\rm ext} =( \mathbf{S}_{i}^{\alpha} -\frac{m_{\alpha}}{M}\mathbf{S} +k_BT\frac{M_{tot}-Nm_{\alpha}}{M_{tot}N}\mathbf{I})\cdot \boldsymbol{F}_{\rm e}\]

where \(\mathbf{S}_{i}^{\alpha}=E_{i}^{\alpha}\mathbf{I}+\mathbf{W}_{i}^{\alpha}\), \(\mathbf{S}=\sum_{\beta=1}^{M}\sum_{i=1}^{N_{\beta}}\mathbf{S}_{i}^{\beta}\), \(M_{tot}\) is the total mass of the system, \(N\) is the atom number of the system. Any physical quantity \(A(t)\) is related to driving force by correlation funtion:

\[\langle \boldsymbol{A}(t) \rangle =\langle A(0) \rangle + (\int_{0}^{t}dt'\frac{\langle \boldsymbol{A}(t') \otimes \boldsymbol{J}_{q}(0) \rangle}{k_BT})\cdot \boldsymbol{F}_{\rm e}\]

In the case of momentum flux \(\boldsymbol{J}_{\gamma}\) of \(\gamma\) component as dissipative flux:

\[\boldsymbol{F}_{i}^{\alpha,\rm ext}=c_{\alpha}\boldsymbol{F}_{\rm e}\]

where \(c_{\gamma}=\frac{N}{N_{\gamma}}\), and \(c_{\beta}=-\frac{Nm_{\beta}}{M'}\), \(M'=\sum_{\epsilon=1,\epsilon\neq\gamma}^{M}N_{\epsilon}m_{\epsilon}\). Similar to the former case,

\[\langle \boldsymbol{A}(t) \rangle =\langle A(0) \rangle + (\frac{N}{M'}+\frac{N}{N_{\gamma}m_{\gamma}})(\int_{0}^{t}dt'\frac{\langle \boldsymbol{A}(t') \otimes \boldsymbol{J}_{\gamma}(0) \rangle}{k_BT})\cdot \boldsymbol{F}_{\rm e}\]

Then we can obtain any matrix element \(\Lambda_{ij}\) of onsager matrix by:

\[\Lambda_{ij} =\frac{1}{k_BV}\int_{0}^{t}dt'\langle \boldsymbol{J}_{i}(t') \otimes \boldsymbol{J}_{j}(0) \rangle\]

The onsager matrix is arranged as:

\[\begin{split}\begin{array}{cccc} \Lambda_{qq} & \Lambda_{q1} & \cdots & \Lambda_{q(M-1)} \\ \Lambda_{1q} & \Lambda_{11} & \cdots & \Lambda_{1(M-1)} \\ \vdots & \vdots & \ddots & \vdots \\ \Lambda_{(M-1)q} & \Lambda_{(M-1)1} & \cdots & \Lambda_{(M-1)(M-1)} \end{array}\end{split}\]

The thermal conductivity could be derived from onsager matrix:

\[\kappa=\frac{1}{T^{2}(\Lambda^{-1})_{00}}\]

The compute_hnemdec keyword should be used to add the driving force and calculate the thermal conductivity.
The computed results are saved to the onsager.out file.