HNEMD Thermal Transport

This example is originally from GPUMD and has been added here with minor changes only to demonstrate how to use gpyumd.

1. Introduction

• In this tutorial, we use the HNEMD method to study heat transport in graphene at 300 K and zero pressure. Here we consider the diffusive regime. The ballistic regime can be found in the NEMD example. The spectral decomposition method as described in [Fan 2019] is used here.

Importing Relevant Functions

• The inputs/outputs for GPUMD are processed using the Atomic Simulation Environment (ASE) and the gpyumd package.

import matplotlib.pyplot as plt
import numpy as np
from gpyumd.atoms import GpumdAtoms
from gpyumd.load import load_shc, load_kappa
from ase.build import graphene_nanoribbon
from gpyumd.sim import Simulation
from gpyumd.math import running_ave
from gpyumd.calc import calc_spectral_kappa
import gpyumd.keyword as kwd

2. Preparing the Inputs

• We consider a graphene sheet of size of about 25 nm x 25 nm (24000 atoms). The transport is in the \(y\) direction. We divide the length in the \(y\) direction into 2 groups, with 4000 atoms in group 0 and 20000 atoms in group 1

• We use a Tersoff potential [Tersoff 1989] parameterized by Lindsay et al. [Lindsay 2010].

• Note that the xyz.in file will be written at the same time as the run.in file with the simulation object below.

gnr = GpumdAtoms(graphene_nanoribbon(100, 60, type='armchair', sheet=True, vacuum=3.35/2))
gnr.euler_rotate(theta=90)
lx, lz, ly = gnr.cell.lengths()
gnr.cell = gnr.cell.new((lx, ly, lz))
gnr.center()
gnr.pbc = [True, True, False]
gnr.set_cutoff(2.1)
gnr.set_max_neighbors(3)
gnr

GpumdAtoms(symbols='C24000', pbc=[True, True, False], cell=[245.95121467478057, 255.6, 3.35])

Add Groups

Ly = 3*1.42*10
split = [0, Ly, 255.6]
group_method, ncounts = gnr.group_by_position(split, direction='y')
print("Atoms per group:", ncounts)
print("Total atoms:", sum(ncounts))

Atoms per group: [ 4000 20000]
Total atoms: 24000

The run.in file

The gpyumd package can be used to generate valid run.in input files as well as other necessary input files. It follows the definitions described in the inputs and outputs documentation for GPUMD.

hnemd_sim = Simulation(gnr, driver_directory='.')

# Equilibration run
header_comment = "Equilibration"
curr_run = hnemd_sim.add_run(number_of_steps=1e6, run_name='equilibration', run_header=header_comment)

keywords = [
    kwd.Velocity(initial_temperature=300),
    kwd.EnsembleNVT(method='nvt_nhc', initial_temperature=300,
                    final_temperature=300, thermostat_coupling=100),
    kwd.NeighborOff(),
    kwd.TimeStep(dt_in_fs=1),
    kwd.DumpThermo(1000)
]

for keyword in keywords:
    curr_run.add_keyword(keyword)

# Production run
header_comment = "Production"
curr_run = hnemd_sim.add_run(run_name='production', run_header=header_comment)

# Time step is propagated from the last run
keywords = [
    kwd.EnsembleNVT(method='nvt_nhc', initial_temperature=300,
                    final_temperature=300, thermostat_coupling=100),
    kwd.NeighborOff(),
    kwd.ComputeHNEMD(output_interval=1000, driving_force_y=1e-5),  # other driving forces default to 0
    kwd.ComputeSHC(sample_interval=2, num_corr_steps=250, transport_direction='y', num_omega=1000, \
                   max_omega=400, grouping_method=0, group_id=0),
    kwd.RunKeyword(number_of_steps=1e7)  # Can enter run keyword here instead of in the 'add_run' method
]

for keyword in keywords:
    curr_run.add_keyword(keyword)

potential_directory = "/path/to/GPUMD/potentials/tersoff"
tersoff_potential = \
    kwd.Potential(filename='Graphene_Lindsay_2010_modified.txt', symbols=['C'], directory=potential_directory)
hnemd_sim.add_potential(tersoff_potential)

hnemd_sim.create_simulation(copy_potentials=True)

The run.in input file is given below:

potential Graphene_Lindsay_2010_modified.txt 0

# Equilibration
velocity 300
time_step 1
ensemble nvt_nhc 300 300 100
neighbor off
dump_thermo 1000
run 1000000

# Production
ensemble nvt_nhc 300 300 100
neighbor off
compute_hnemd 1000 0 1e-05 0
compute_shc 2 250 1 1000 400 group 0 0
run 10000000

• The first line uses the potential keyword to define the potential to be used, which is specified in the file Graphene_Lindsay_2010_modified.txt.

• The second line uses the velocity keyword and sets the velocities to 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 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 (ensemble) with the same parameters as in the equilibration stage.

• The compute_hnemd is used to add a driving force and compute the thermal conductivity using the HNEMD method [Fan 2019]. Here, the conductivity data will be averaged for each 1000 steps before written out, and the driving force parameter is \(10^{-5}\) A-1 and is in the \(y\) direction.

• The line with the compute_shc keyword is used to compute the spectral heat current (SHC). The SHC will be calculated for group 0 in 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). The number of frequency points is 1000 and the maximum angular frequncy is 400 THz.

• There are \(10^7\) steps (10 ns) in the production run. This is just an example. To get more accurate results, we suggest you use \(2 \times 10^7\) steps (20 ns) and do a few independent runs and then average the relevant data.

3. Results and Discussion

Computation Time

• The HNEMD simulation takes about 2.5 hours using a laptop with a GeForce RTX 2070 GPU and about 38 minutes with a Tesla V100.

aw = 2
fs = 16
font = {'size'   : fs}
plt.rc('font', **font)
plt.rc('axes' , linewidth=aw)

def set_fig_properties(ax_list):
    ax_list = ax_list if isinstance(ax_list, list) else [ax_list]
    tl = 8
    tw = 2
    tlm = 4

    for ax in ax_list:
        ax.tick_params(which='major', length=tl, width=tw)
        ax.tick_params(which='minor', length=tlm, width=tw)
        ax.tick_params(which='both', axis='both', direction='in', right=True, top=True)

Plot HNEMD Results

• The kappa.out output file is loaded and processed to create the following figure.

kappa = load_kappa()
kappa.keys()

dict_keys(['kxi', 'kxo', 'kyi', 'kyo', 'kz'])

t = np.arange(1,kappa['kxi'].shape[0]+1)*0.001  # ns
kappa['kyi_ra'] = running_ave(kappa['kyi'],t)
kappa['kyo_ra'] = running_ave(kappa['kyo'],t)
kappa['kxi_ra'] = running_ave(kappa['kxi'],t)
kappa['kxo_ra'] = running_ave(kappa['kxo'],t)
kappa['kz_ra'] = running_ave(kappa['kz'],t)

axes = list()

plt.figure(figsize=(12,10))
plt.subplot(2,2,1)
axes.append(plt.gca())
plt.plot(t, kappa['kyi'],color='C7',alpha=0.5)
plt.plot(t, kappa['kyi_ra'], linewidth=2)
plt.xlim([0, 10])
plt.gca().set_xticks(range(0,11,2))
plt.ylim([-2000, 4000])
plt.gca().set_yticks(range(-2000,4001,1000))
plt.xlabel('time (ns)')
plt.ylabel(r'$\kappa_{in}$ W/m/K')
plt.title('(a)')

plt.subplot(2,2,2)
axes.append(plt.gca())
plt.plot(t, kappa['kyo'],color='C7',alpha=0.5)
plt.plot(t, kappa['kyo_ra'], linewidth=2, color='C3')
plt.xlim([0, 10])
plt.gca().set_xticks(range(0,11,2))
plt.ylim([0, 4000])
plt.gca().set_yticks(range(0,4001,1000))
plt.xlabel('time (ns)')
plt.ylabel(r'$\kappa_{out}$ (W/m/K)')
plt.title('(b)')

plt.subplot(2,2,3)
axes.append(plt.gca())
plt.plot(t, kappa['kyi_ra'], linewidth=2)
plt.plot(t, kappa['kyo_ra'], linewidth=2, color='C3')
plt.plot(t, kappa['kyi_ra']+kappa['kyo_ra'], linewidth=2, color='k')
plt.xlim([0, 10])
plt.gca().set_xticks(range(0,11,2))
plt.ylim([0, 4000])
plt.gca().set_yticks(range(0,4001,1000))
plt.xlabel('time (ns)')
plt.ylabel(r'$\kappa$ (W/m/K)')
plt.legend(['in', 'out', 'total'])
plt.title('(c)')


plt.subplot(2,2,4)
axes.append(plt.gca())
plt.plot(t, kappa['kyi_ra']+kappa['kyo_ra'],color='k', linewidth=2)
plt.plot(t, kappa['kxi_ra']+kappa['kxo_ra'], color='C0', linewidth=2)
plt.plot(t, kappa['kz_ra'], color='C3', linewidth=2)
plt.xlim([0, 10])
plt.gca().set_xticks(range(0,11,2))
plt.ylim([0, 4000])
plt.gca().set_yticks(range(-2000,4001,1000))
plt.xlabel('time (ns)')
plt.ylabel(r'$\kappa$ (W/m/K)')
plt.legend(['yy', 'xy', 'zy'])
plt.title('(d)')

set_fig_properties(axes)
plt.tight_layout()
plt.show()

(a) In-plane thermal conductivity \(\kappa_{\rm in}(t)\) as a function of production time. (b) Out-of-plane thermal conductivity \(\kappa_{\rm out}(t)\) as a function of production time. (c) In-plane, out-of-plane, and total thermal conductivity as a function of production time. (d) Comparing \(\kappa_{yy}\), \(\kappa_{xy}\), and \(\kappa_{zy}\).

• 1. The gray line represents the instant in-plane thermal conductivity \(\kappa_{\rm in}(t)\) from the kappa.out file. The dashed line represents the cumulative time average:

     \[\kappa^{\rm ave}_{\rm in}(t) = \frac{1}{t} \int_0^t \kappa_{\rm in}(t') dt'\]

• 2. Similar to (a), but for the out-of-plane thermal conductivity \(\kappa_{\rm out}(t)\) and \(\kappa^{\rm ave}_{\rm out}(t)\).

• 3. Cumulative time averaged in-plane, out-of-plane, and total thermal conductivity: \(\kappa^{\rm ave}_{\rm in}(t)\), \(\kappa^{\rm ave}_{\rm out}(t)\) and \(\kappa^{\rm ave}_{\rm total}(t) = \kappa^{\rm ave}_{\rm in}(t) + \kappa^{\rm ave}_{\rm out}(t)\). It is clear that the out-of-plane phonons have a much large average phonon relaxation time (scattering time) than the in-plane phonons, and the thermal conductivity in pristine graphene is dominated by out-of-plane (flexural) phonons. This is consistent with the results from the EMD simulation as discussed in the EMD tutorial.

• 4. Some cumulative time averaged thermal conductivity components of the tensor: \(\kappa^{\rm ave}_{yy}(t)\), \(\kappa^{\rm ave}_{xy}(t)\) and \(\kappa^{\rm ave}_{zy}(t)\). It is clear that the off-diagonal components of the thermal conductivity tensor are zero.

Plot Spectral Heat Current Results

• The shc.out output file is loaded and processed to create the following figure.

shc = load_shc(num_corr_points=250, num_omega=1000)['run0']
shc.keys()

dict_keys(['t', 'Ki', 'Ko', 'nu', 'jwi', 'jwo'])

l = gnr.cell.lengths()
Lx, Lz = l[0], l[2]
V = Lx * Ly * Lz
T = 300
Fe = 1.0e-5
calc_spectral_kappa(shc, driving_force=Fe, temperature=T, volume=V)
shc['kw'] = shc['kwi'] + shc['kwo']
shc['K'] = shc['Ki'] + shc['Ko']
Gc = np.load('Gc.npy')
shc.keys()

dict_keys(['t', 'Ki', 'Ko', 'nu', 'jwi', 'jwo', 'kwi', 'kwo', 'kw', 'K'])

lambda_i = shc['kw']/Gc
length = np.logspace(1,6,100)
k_L = np.zeros_like(length)
for i, el in enumerate(length):
    k_L[i] = np.trapz(shc['kw']/(1+lambda_i/el), shc['nu'])

axes = list()

plt.figure(figsize=(12,10))
plt.subplot(2,2,1)
axes.append(plt.gca())
plt.plot(shc['t'], shc['K']/Ly, linewidth=3)
plt.xlim([-0.5, 0.5])
plt.gca().set_xticks([-0.5, 0, 0.5])
plt.ylim([-1, 5])
plt.gca().set_yticks(range(-1,6,1))
plt.ylabel('K (eV/ps)')
plt.xlabel('Correlation time (ps)')
plt.title('(a)')

plt.subplot(2,2,2)
axes.append(plt.gca())
plt.plot(shc['nu'], shc['kw'],linewidth=3)
plt.xlim([0, 50])
plt.gca().set_xticks(range(0,51,10))
plt.ylim([0, 200])
plt.gca().set_yticks(range(0,201,50))
plt.ylabel(r'$\kappa$($\omega$) (W/m/K/THz)')
plt.xlabel(r'$\nu$ (THz)')
plt.title('(b)')

plt.subplot(2,2,3)
axes.append(plt.gca())
plt.plot(shc['nu'], lambda_i,linewidth=3)
plt.xlim([0, 50])
plt.gca().set_xticks(range(0,51,10))
plt.ylim([0, 6000])
plt.gca().set_yticks(range(0,6001,1000))
plt.ylabel(r'$\lambda$($\omega$) (nm)')
plt.xlabel(r'$\nu$ (THz)')
plt.title('(c)')

plt.subplot(2,2,4)
axes.append(plt.gca())
plt.semilogx(length/1000, k_L,linewidth=3)
plt.xlim([1e-2, 1e3])
plt.ylim([0, 3000])
plt.gca().set_yticks(range(0,3001,500))
plt.ylabel(r'$\kappa$ (W/m/K)')
plt.xlabel(r'L ($\mu$m)')
plt.title('(d)')

set_fig_properties(axes)
plt.tight_layout()
plt.show()

(a) Virial-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 [Fan 2019].

• 1. The virial-velocity correlation function \(K(t)\). See Theoretical formulations for the definition of this quantity.

• 2. The spectral thermal conductivity \(\kappa(\omega)\). See Theoretical formulations for the definition of this quantity.

• 3. The spectral phonon mean free path calculated as [Fan 2019]

     \[ \begin{align}\begin{aligned} \lambda(\omega)=\kappa(\omega)/G(\omega),\\where :math:`G(\omega)` is the quasi-ballistic spectral thermal conductance from the above NEMD simulation.\end{aligned}\end{align} \]

• 4. The length-dependent thermal conductivity calculated as [Fan 2019]

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

4. References

• [Fan 2019] 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).

• [Lindsay 2010] L. Lindsay and D.A. Broido, Optimized Tersoff and Brenner emperical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene, Phys. Rev. B, 81, 205441 (2010).

• [Tersoff 1989] J. Tersoff, Modeling solid-state chemistry: Interatomic potentials for multicomponent systems, Phys. Rev. B 39, 5566(R) (1989).