Use HNEMD Method to Calculate Thermal Conductivity

Calculating Thermal Conductivity with the HNEMD Method

Method Overview

In materials science and engineering, thermal conductivity is a key parameter for measuring a material's ability to conduct heat, and it directly affects heat dissipation in electronic devices, the design of energy materials, and the application of nanotechnology. Traditional experimental measurement methods are often time-consuming and costly, whereas molecular dynamics (MD) simulation is a powerful theoretical tool that can reveal heat transport mechanisms at the atomic scale. However, traditional non-equilibrium molecular dynamics (NEMD) and equilibrium molecular dynamics (EMD) methods often suffer from boundary effects and low computational efficiency in simulations. The HNEMD (Homogeneous Non-Equilibrium Molecular Dynamics) method applies a homogeneous heat flux rather than relying on a temperature gradient, thereby avoiding the boundary perturbation caused by temperature gradients in traditional methods. It is particularly suitable for nanomaterials, low-dimensional systems, and other boundary-sensitive scenarios, significantly improving the efficiency and accuracy of thermal conductivity calculations. Moreover, because the system temperature is uniform, it is, like EMD, insensitive to system size. Compared with the EMD method, the HNEMD method improves computational efficiency by two orders of magnitude.

Algorithm Overview

Below is an explicit algorithm, proposed in the literature, that can be implemented directly on a computer. An HNEMD simulation consists of the following steps:

1. Equilibration stage. First, as in any molecular dynamics (MD) simulation, the system is equilibrated in the $NVT$ or $NPT$ ensemble to reach thermal equilibrium. Note that, as with the EMD method, periodic boundary conditions must be applied along the transport direction.

2. Simulation stage. Second, a uniform heat flux is generated by superimposing on the interatomic forces the driving force given by

$${\mathbf{D}}_i\cdot {\mathbf{F}}_e=E_i{\mathbf{F}}_e+\sum _{j\ne i}(\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}\otimes {\mathbf{r}}_{ij})\cdot {\mathbf{F}}_e.$$$${\mathbf{F}}_i=\sum _{j\ne i}(\frac{\partial U_i}{\partial {\mathbf{r}}_{ij}}-\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}),$$

which gives the total force:

$${\mathbf{F}}_i^{\text{tot}}={\mathbf{F}}_i+{\mathbf{D}}_i\cdot {\mathbf{F}}_e.$$

To conserve the total momentum of the system, the average force must be subtracted from the force on each particle.

$${\mathbf{F}}_i^{\text{tot}}\to {\mathbf{F}}_i^{\text{tot}}-\frac{1}{N}\sum _i{\mathbf{F}}_i^{\text{tot}}.$$

At this stage, a thermostat must also be applied to keep the system temperature at the target value; otherwise, the system would heat up because of the driving force.

3. Post-processing. Finally, we sample the heat flux according to

$${\mathbf{J}}_q={\mathbf{J}}_q^{\text{kin}}+{\mathbf{J}}_q^{\text{pot}}=\sum _i\frac{{\mathbf{p}}_i}{m_i}{E}_i+\sum _{i,j\ne i}\frac{{\mathbf{p}}_i}{m_i}\cdot (\frac{\partial U_j}{\partial {\mathbf{r}}_{ji}}\otimes {\mathbf{r}}_{ij})$$

and compute the thermal conductivity from

$$\kappa (t)=\frac{1}{t}{\int }_0^{t}ds\frac{⟨J_q(s){⟩}_{\text{ne}}}{TV{F}_e}$$

LAMMPS Implementation

After several rounds of trial and error, I implemented a number of different approaches, and finally managed to carry out the entire calculation within a LAMMPS script. Of course, a LAMMPS plugin version also exists, though I do not know whether anyone would be interested. It took a great deal of effort, and in the end it was not of much use......

这里好像有些什么

The example below is a LAMMPS input file that computes the thermal conductivity of the KA system.

# =======================================================================
# ================================ System ===============================
units lj
atom_style atomic
dimension 3
boundary p p p
# ============================== set System =============================
variable    Temp        equal   0.8
variable    dt          equal   0.005
variable    natoms      equal   1000
variable    density     equal   1.2
variable    volume      equal   ${natoms}/${density}
variable    box_length  equal   ${volume}^(1.0/3.0)
# ============================= restart part ============================
variable    restart_file_exists equal 1
restart 1000000 restart.*.restart
if "${restart_file_exists} == 1" then &
  "read_restart restart.1000000.restart" &
  "jump SELF HNEMD"
# ============================== create box =============================
region box block 0 ${box_length} 0 ${box_length} 0 ${box_length}
create_box 2 box
# ============================= create atoms ============================
variable    natoms1     equal   round(${natoms}*0.8)
variable    natoms2     equal   ${natoms}-${natoms1}
create_atoms 1 random ${natoms1} 12345 box
create_atoms 2 random ${natoms2} 67890 box
mass 1 1.0
mass 2 1.0
# ========================== create Potential ===========================
# Lennard-Jones Potential paramaters
pair_style  lj/cut  2.5
pair_coeff  1   1   1.0 1.0     2.5
pair_coeff  1   2   1.5 0.8     2.0
pair_coeff  2   2   0.5 0.88    2.2
# ============================ create group =============================
group A type 1
group B type 2
# ============================== MD parts ===============================
# ============================ neighbor list ============================
neighbor 0.5 bin
neigh_modify delay 0 every 3 check yes
# ============================== minimize ===============================
min_style sd
minimize 1.e-7 1.e-7 1000 5000
# write_data "minnimize.data" pair ij
# ============================== initilize ==============================
velocity all create ${Temp} 188546 dist gaussian mom yes loop all
timestep ${dt}
# ============================ equilibrium ==============================
fix 1 all nvt temp ${Temp} ${Temp} 0.5
thermo 100
thermo_style custom step temp pe etotal
run 1000000
# =======================================================================
label HNEMD
# thermal conductivity calculation
# =======================================================================
# ======== Homogeneous nonequilibrium molecular dynamics method =========
# ============================== Parameters =============================
variable    Fe_ix   equal   0.005
variable    Fe_iy   equal   0.005
variable    Fe_iz   equal   0.005
variable    Fe_x    atom    v_Fe_ix
variable    Fe_y    atom    v_Fe_iy
variable    Fe_z    atom    v_Fe_iz
# =======================================================================
# =============================== Scripts ===============================
compute     myKE    all     ke/atom
compute     myPE    all     pe/atom
compute     Wij     all     stress/atom NULL virial
# ===============================  begin  ===============================
# ====================== compute and add the force ======================
variable    Eng     atom    c_myKE+c_myPE
variable    fex     atom    v_Eng*v_Fe_x-c_Wij[1]*v_Fe_x-c_Wij[4]*v_Fe_y-c_Wij[5]*v_Fe_z
variable    fey     atom    v_Eng*v_Fe_y-c_Wij[4]*v_Fe_x-c_Wij[2]*v_Fe_y-c_Wij[6]*v_Fe_z
variable    fez     atom    v_Eng*v_Fe_x-c_Wij[5]*v_Fe_x-c_Wij[6]*v_Fe_y-c_Wij[3]*v_Fe_z
variable    Jex     atom    v_Eng*vx-c_Wij[1]*vx-c_Wij[4]*vy-c_Wij[5]*vz      # debug
variable    Jey     atom    v_Eng*vy-c_Wij[4]*vx-c_Wij[2]*vy-c_Wij[6]*vz      # debug
variable    Jez     atom    v_Eng*vz-c_Wij[5]*vx-c_Wij[6]*vy-c_Wij[3]*vz      # debug
compute     hc  all reduce  sum v_Jex   v_Jey   v_Jez
compute     Fe  all reduce  sum v_fex   v_fey   v_fez
# ===================== reduce force of mass center =====================
variable    fxe     equal  -c_Fe[1]
variable    fye     equal  -c_Fe[2]
variable    fze     equal  -c_Fe[3]
variable    fix_fx  atom    v_fxe/${natoms}
variable    fix_fy  atom    v_fye/${natoms}
variable    fix_fz  atom    v_fze/${natoms}
variable    add_fx  atom    v_fix_fx+v_fex
variable    add_fy  atom    v_fix_fy+v_fey
variable    add_fz  atom    v_fix_fz+v_fez
fix Fe  all addforce    v_add_fx    v_add_fy    v_add_fz
# ============================  heat flux  ==============================
compute     flux    all     heat/flux   myKE    myPE    Wij    
variable    Jx      equal   c_flux[1]
variable    Jy      equal   c_flux[2]
variable    Jz      equal   c_flux[3]
variable    Jhx     equal   c_flux[4]
variable    Jhy     equal   c_flux[5]
variable    Jhz     equal   c_flux[6]
# ========================  Thermal Conductivity  =======================
variable    K_x     equal   v_Jx/(${Temp}*vol*v_Fe_ix)
variable    K_y     equal   v_Jy/(${Temp}*vol*v_Fe_iy)
variable    K_z     equal   v_Jz/(${Temp}*vol*v_Fe_iz)
variable    K_sum   equal   (v_K_x+v_K_y+v_K_z)/3
# =====================  Thermal Conductivity dump  =====================
fix kappa   all print   1   "${K_x} ${K_y} ${K_z} ${K_sum}" file    data.kappa  screen  no 
# ==========================  heat flux dump  ===========================
# dump heatflux all custom 1000 dump.custom id type v_Jex v_Jey v_Jez
### dump  debug all custom 100 dump.velocities id type fx fy fz v_fex v_fey v_fez v_Eng  # debug
# ======================  end of the modify part  =======================
# =======================================================================
# =======================================================================
# ================================== MD =================================
fix 1 all nvt temp ${Temp} ${Temp} 0.5
thermo 100
thermo_style custom step temp pe etotal v_K_x v_K_y v_K_z v_K_sum v_Jx v_Jy v_Jz
# thermo_modify norm no
run 1000000