Use HNEMD Method to Calculate Thermal Conductivity

使用HNEMD方法计算热导率

方法简介

在材料科学和工程领域，热导率是衡量材料导热性能的关键参数，直接影响电子器件散热、能源材料设计和纳米技术应用。传统的实验测量方法往往耗时且成本高昂，而分子动力学（MD）模拟作为一种强大的理论工具，能够从原子尺度揭示热输运机制。然而，传统的非平衡分子动力学（NEMD）方法与平衡分子动力学（EMD）方法在模拟中常面临边界效应和计算效率低的问题。HNEMD（Homogeneous Non-Equilibrium Molecular Dynamics） 方法通过均匀施加热流而非依赖温度梯度，避免了传统方法中因温度梯度导致的边界扰动问题，尤其适用于纳米材料、低维体系等对边界敏感的场景，显著提升了热导率计算的效率和准确性，同时由于体系温度均一，从而与EMD一样对体系大小不敏感。HNEMD方法相比于EMD方法在计算效率上提升了两个数量级。

算法简介

根据文献提出的一种可以直接在计算机上实现的明确算法。HNEMD模拟包括以下步骤：

1. 平衡阶段。首先，与任何分子动力学（MD）模拟一样，我们在 或 系综中对系统进行平衡，以达到热平衡。需要注意的是，与EMD方法相同，必须在输运方向上应用周期性边界条件。

2. 模拟阶段。其次，通过在原子间作用力上叠加由公式

给出的驱动力来生成均匀的热流：

从而得到总力：

为了保持系统的总动量守恒，必须从每个粒子的力中减去系统的平均力。

在此阶段，还需要应用热浴以保持系统温度在目标值；否则，系统会因驱动力而升温。

3. 后处理。最后，我们根据公式

对热流进行采样，并根据公式

计算热导率。

lammps实现

经过几轮折腾，实现过多种方法，最后实现了在LAMMPS脚本中实现全部计算。 当然lammps插件形式也有，不知道有没有人感兴趣。 花了很多力气高，最后没有派上什么用场......

下面的例子实现了计算KA体系热导率的lammps输入文件。

# =======================================================================
# ================================ 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