随着现代电子器件尺寸的不断下降，我们可以看到纳米器件的复杂性会不断增加。由两种以上材料组成的器件的热导率是很难预测的。在大多数情况下，器件工程师想要利用不同材料构成界面来增大或减小热导率。

在本教程中，你将学习如何利用非平衡态分子动力学结合经验力场（势）的方法模拟热流通过界面，计算界面处的热导（热阻）。

通过在热源和冷源间施加温度梯度，可以模拟经硅晶体的晶界的热流。

在给定温度梯度下，块体材料的热导率 $K_{bulk}$ 与热流的关系由傅立叶方程给出：

$\frac{\langle dQ/dt\rangle}{A} = -\kappa_{bulk}\langle dT/dz\rangle$

在这里， ⟨dQ/dt⟩ 代表沿着温度梯度的平均热流（也就是单位时间传输的动能），A 是垂直于热流方向的横截面积。

处在热源和冷源中间的界面通常会使温度分布产生一个 ΔT 的突变（温度的“跳跃”），因为界面会对热流产生额外的热阻 R。在这种情况下，界面处热导 G，通常被称作卡皮查（Kapitza）热导，被定义为

$G=\frac{1}{R}=\frac{\langle dQ/dt\rangle}{A\Delta T}$

界面导热性质可以通过分子动力学模拟进行计算。最普遍的方法是在分子动力学模拟中引入温度梯度，也就是所谓的非平衡动量交换或反向非平衡分子动力学方法¹⁾²⁾。在这种方法下，两端的两组原子被分别选择作为热源和冷源。在模拟期间，冷源中最热的原子和热源中最冷的原子根据特殊的标准进行动量交换，保证总能量和动量守恒³⁾。这种外部的热流则通过系统内由热源流向冷源的热流得到平衡，使体系内产生可观测的温度梯度（以及有界面时产生的温度的跳跃）。

这里考查一个由硅单晶的两种不同晶向形成的晶界模型。在模拟中，虽然你可以让热流沿着任意方向。但是只有沿着z方向的温度梯度可以求算。因此，构建体系时必须确保界面结构指向正确的方向。必要时可以通过旋转体系确保热流方向沿着z轴。

下面在硅的（100）表面和（110）表面之间创建界面。打开 Builder，通过 Add ‣ From database 添加 alpha-silicon，并且：

• 生成（100）表面基胞，打开 Bulk Tools ‣ Supercell plugin，点击 Conventional and Transform
• 生成（110）表面基胞
  • 打开位于右手边 Builders 下的 Surface (Cleave) 工具
  • 选择（110）米勒指数，点击 Next，保留 1×1 面晶格，点击 Next，选择 Periodic (bulk-like)
  • 点击 Finish 完成表面创建

• 生成（100）表面基胞，打开 Bulk Tools ‣ Supercell plugin, 点击 Conventional and Transform.
• 在（110）结构中重复上面的步骤

为了避免因连接两种不同晶格常量的单胞产生太大的形变，需要选择合适的重复以获得相称的单胞尺寸。单胞尺寸可以在 Bulk tools ‣ Lattice parameters 中观察并修改。

（100）面单胞的尺寸为 5.43 x 5.43 x 5.43 $Å^3$, ，（110）面 为3.84 x 5.43 x 3.84 $Å^3$。因此，（100）面的单胞应该沿着 ABC 分别进行 5x5x25 次的重复，（110）则需要 7x5x35 的重复，以使提供的超胞有着几乎相等的尺寸。你可以用 Bulk tools ‣ Repeat 插件去调节重复。设置想要的重复次数并点击Apply.

创建一个界面可能是一件繁琐而枯燥的事情，这时可以使用 Interface Builder 帮你更简单灵活的完成这个建模过程。打开右边的 Builders > Interface 选项，将（100）超胞拖拽到第一个方框中，将 (110) 超胞放到第二个方框中。

界面的预览图应该会自动出现（如果没有则点击 Preview）。你将注意到在两个表面间有一条缝隙。你可以通过点击 Shift Surfaces 合拢该缝隙，在 z 方向填入一个负的偏移量直到两个表面间出现化学键。不仅如此，你也可以沿 y 方向调整对齐两表面。自动生成的默认值作为初始结构一般是不错的。

在此计算中超胞在 c 方向上是非周期的，这意味着系统上下两端被两表面截断，热流仅能沿着 z 轴的正方向。横向尺寸（A和B）保留周期性。

为获得这样的设置，可以简单地增加 c 方向的长度。你可以通过打开 Bulk Tools ‣ Lattice Parameters，将 Cz 项由 268 Å 增加到 300 Å 来实现。你需要确保用的是笛卡尔坐标。

最后，片层结构应该置于超胞的中间。因此，打开 Coordinate Tools ‣ Center，不勾选 A 和 B，点击 Apply.

在模拟实际热流之前，体系结构必须要达到很好的平衡结构，以减小因构造界面人为产生的缺陷，去除横向应力；最后得到一个在目标温度下平衡的体系。为了达到这个目的，你需要通过点击 Send 图标将最终确定的界面发送到Script Generator ，并添加 3 个模块:

• New Calculator
• Optimize ‣ OptimizeGeometry
• Optimize ‣ MolecularDynamics

将默认的输出文件名改为 Silicon_grain_boundary.nc，并打开 New Calculator 模块

• 选择 ATK-Classical 和 StillingerWeber_Si_1985 势。
• 不勾选 save 和 print，通过点击 OK 关闭小窗口。

打开 OptimizeGeometry 模块，将最大力参数调到 0.1 eV/Å 作为阈值。其余参数保留原默认值，关闭小窗口。

接下来的步骤是用分子动力学手段让横向晶胞矢量得到弛豫、平衡，使之达到目标温度。整个热输运模拟过程要在平均 600K 的温度下进行。为达到这个目的，打开 MolecularDynamics 模块，将参数设置成 NPT 模拟，这是一种在常温常压下的模拟。

• 将类型选为 NPT Melchionna。这种方法将各个晶胞矢量看成独立的，这对和本例一样的各向异性体系来说是非常有必要的（你也可以看看分子动力学基础）。
• 设置 number of steps 为 50000 步，log interval 为 500.
• 选择一个合适的 trajectory 名称，比如， GrainBoundary_equilibrate_NPT_traj.nc
• 将热库设置成 600K，外部应力为 0 bar。
• 将恒温时间刻度和恒压时长设置成 50.0 fs，并关闭 zz 方向的应力耦合，因为这个方向的自由表面会自动被调整。
• 选择 Maxwell-Boltzmann 分布作为初始速度，初始温度为 600K，点击方框移除质心动量。
• 不勾选 save 和 print，因为 trajectory 被保存了

最后设置应该如下所示：

将脚本发送到 Editor 来增加一些额外的功能，因为你必须在结果中提取 A- 和 B-向量的平均长度。拖到底部，在脚本的最底下加上下列几行。

# Get the trajectory of cell vectors.
cells = md_trajectory.cells()
# Extract the xx- and yy-components.
vector_xx = [cell[0, 0].inUnitsOf(Angstrom) for cell in cells]
vector_yy = [cell[1, 1].inUnitsOf(Angstrom) for cell in cells]
# Calculate and print the mean value.
print "Average xx-component:  ", numpy.mean(vector_xx)
print "Average yy-component:  ", numpy.mean(vector_yy)

将脚本发送到 JobManager 并运行。

在模拟结束后，A 和 B 的平均晶格矢量将会在 log 文件中的底端被打印出来。选择 Lab floor 中的 trajectory 并在 Movie Tool 中打开它。播放运动轨迹并将其最后一帧通过 标志发送到 Builder 中。

在 Builder，打开 Bulk Tools ‣ Lattice parameters 对话框。选择 Keep Fractional ，将晶胞矢量中的 xx- 和 yy-元修改成从模拟中获得的平均值。关闭 Lattice parameters对话框。在以后的模拟中，晶格常数将会被固定。

现在你必须定义热流模拟过程中的热源和冷源区域，放大结构最左端，选择最外面的8个原子层（即两个单胞），如下图所示。

• 打开右侧工具 Select ‣ Tags ，为标签（tag）输入一个合适的名字，例如，heat_source，并按回车.被选中的原子将被标记为这个标签名字，并且能在模拟过程能被引用。
• 在结构的另一端重复相同的过程。这次用一个不同的标签（tag）名称，比如，heat_sink。
• 将该结构发送到 Script Generator ，添加 New Calculator 和两个分子动力学模块。
• 按上述同样的方法编辑计算器模块，即：选择 ATK-Classical 和 StillingerWeber_Si_1985 势。

首先，你运行另一个处在 600K 温度下的平衡模拟，这次在恒定体积下准备一个平衡态体系。打开第一个分子动力学模块 并，

• 选择类型为 NVT Nose Hoover Chain
• 设置步数为 100000，log interval 为 5000.
• 选择合适的 trajectory 文件，例如，GrainBoundary_eq_NVT_traj.nc
• 选择 600K 温度下 Maxwell-Boltzmann 初始速度，确保已经勾上质心移除选项。
• 设置热库温度为 600K，最终温度也将自动变化到 600K。保留恒温器时间标度 100fs。
• 点击 OK，关闭对话框。

在第二个 MolecularDynamics 模块中，最终的非平衡模拟将被设置。为了实现这一目的：

• 选择类型为 Non-Equilibrium Momentum Exchange
• 设置步数为 800000，长间隔为 1000.
• 选择合适的 trajectory 的文件名，例如，GrainBoundary_nonEq_traj.nc
• 设置初始结构速度，这意味着储存在结构中的速度将被利用。
• 保留时间步长的缺省值为 1fs。
• 将交换间隔项设为 200，这意味着交换每 200 步被触发一次。

你可以使用 tag 里的 “heat_source” 和 “heat_sink” 确定各自热源和冷源区域，这是之前在构建模型时已经定义好的。最终的设置应该如下所示。

• 通过 JobManager 或者保存脚本，在终端中使用 atkpython 运行，对于 8 核集群需要 10 个小时才能完成。

模拟完成后，轨迹文件将出现在 LabFloor 里。用 Movie Tool 打开非平衡轨迹文件，你会看到在整个模拟过程中平均温度是一个常数，片层没有漂移，因为求算总体的平均温度不可能体现温度的空间分布。

用 MD Analyzer 插件打开同一个的轨迹文件，研究热导率，

• 选择温度分布按钮并点击 Load，这可能会耗费一些时间
• 选择计算分布的时间间隔

因此，例如你可以比较 100 ps 与 400 ps 间的分布和 400 ps 和 800 ps 间的分布。用在工具栏中的 Zoom Tool 去放大观察热导分布图中最显著的特征，或者点击 Home 返回初始视图。除了统计性的起伏外，分布应该是一样的。

温度分布显示了一些典型的特征。首先，在分布的首段和末段，可以看到一个正向和一个负向的峰值。这两个部分分别是热源和冷源。这种非线性部分来源于人为的动能转移，这是由动量交换产生的。你应该在你分析的时候去掉这一部分。

第二个显著的特征是在 z=150Å 处，温度会突然跃迁。这是界面区域，因此有额外的热阻。让我们放大该区域更进一步观察一下。

在 T1=625K 和 T2=575K 的间温度变化为 ΔT=50 K。为了计算 Kapitza 热导你还需要知道平均热流。该值可以从日志文件中的 nlprint output 中找到。在该例中，可以获得的热流值 ⟨dQ/dt⟩=0.00191eV/fs。最后你需要的是横截面积，这可以通过 A 和 B 的晶格矢量长度计算得到。根据上面的公式，结合这些数值，就可以计算 Kapitza 热导了。

鉴于单位必须要正确转换，计算热导值最简单的方法是利用 QuantumATK 的内建单位引擎。

• 打开 Editor 并输入下列方程。

dq_dt = 0.00191*eV/fs
A = 26.988*27.118*Angstrom**2
delta_T = 50.0*Kelvin
G = dq_dt / delta_T / A
print "Kapitza conductance in Watt/Kelvin/Meter**2:  ",G.inUnitsOf(Joule/Second/Kelvin/Meter**2)

• 通过 JobManager （或者atkpython）运行脚本，日志文件将显示计算的Kapitza 热导，单位为 Watt/Kelvin/$m^2$。在该例中，结果是 $8.4\times10^8 W/K/m^2$

在该例中，你也可以调节热源或者冷源到边界的距离，大约为 100 Å，来研究晶界尺寸的影响。此时，你需要考虑的是热导率而不是界面热导。用 MD Analyzer 提取温度梯度，在 Zmin 和 Zmax 处输入合理的数值来选择分布的线性区域。

程序将在参数窗口中自动拟合斜率，显示温度梯度。在该例中，（100）取向晶粒的值为 -0.167 K/Å。你对右边的晶粒重复该过程，获得一条相似的斜率。左右两边的斜率不可能一样，这是由于统计起伏和不同的晶粒取向造成的。

现在我们可以在脚本中将计算得到的斜率代入到傅里叶方程中，利用单位时间平均转移的能量、横截面积计算晶粒内的热导率。

dq_dt = 0.00191*eV/fs
A = 26.988*27.118*Angstrom**2
dT_dz = 0.167*Kelvin/Angstrom
G = dq_dt / dT_dz / A
print "Grain conductivity in Watt/Kelvin/Meter:  ",G.inUnitsOf(Joule/Second/Kelvin/Meter)

结果是 25 W/K/m。对于（110）取向晶粒，热导率为 22.5 W/K/m，这些值与模拟的体块状Si的热导率相比是非常的小，在 600K 下，$\kappa_{bulk} = 150 W/K/m$，这是³⁾ 中利用相同的 Stillinger-Weber 势的计算结果。这种结果是可以预料的，因为尺寸受限的晶粒的热导率要低于块体材料中相应的值 ¹⁾ 。这是由于声子的平均自由程通常比 1 μm 要大得多，这意味着模拟的晶粒尺寸抑制了热导^1),3) 。

为了研究晶粒尺寸相关性，你可以重复模拟，准备两种晶胞不同长度的结构，正如参考文献¹⁾ 所示。

块体热导率可以通过多次模拟获得，随着晶粒尺寸l的增加，绘制热导率倒数 1/κ(l) 和晶粒尺寸倒数1/l的关系，并由此外推到无限大晶粒尺寸。

• 英文链接：http://docs.quantumwise.com/tutorials/interfacial_thermal_conductance.html
• 中文翻译：王子群