热导率计算公式：

• k = W/(S⋅Grad(T)) = (dE/dt)/(SΔT/L)

参考文献：

• E.M. Moscarello, el, ACS Appl. Nano Mater. 2022, 5, 10, 13787–13793

模型：

• 在体系中划分一个较高恒温区域为“热源”，类似划分一个较低恒温区域区域为“热漏”，由于二者之间的温差，而产生热传导。由于周期性的缘故，会朝正、反两个方向传递能量，即“热漏”累计接受能量将是两个方向过来的，为了更简单地统计传递能量的速率（即功率 W = dE/dt），需要两个热传导区域大小一样，如此则功率的一半，即为单个热传导区域的传输功率。
• “热漏”接受能量的功率：在AMSinput中，为“热源”、“热漏”分别设置较高、较低温度的恒温系综，两个系综分别排序为系综1，系综2，则计算完毕后，Movie → MD Properties → NHCT1statEnergy、NHCT2statEnergy 代表两个区域的能量随时间而积累的量，那么 NHCT2statEnergy 对时间的梯度的一半，即热传导功率。
• 温度的梯度：温差ΔT、热传导区域的长度L都知道，因此温度梯度 Grad(T) = ΔT/L。
• S指传热区域的横截面面积，在AMSinput → Lattice底部可以直接看到，对于正交Cell来说，如果传热方向为C方向，则S = a×b，以此类推。

确定好热传导区域，剩余的两头就是热源和热漏，如下图所示，红色区域为热源，绿色区域为热漏，热量通过中间环，分别向上、向下传递（注意周期性，如果不明白，请参阅：周期边界条件）

设置普通的分子动力学参数（这里使用的是 UFF 力场，使用 ReaxFF 某个力场或者机器学习力场，类似设置即可）：

注意：

• 步数：设置较长，这里作为演示设置了 20 万步，步数需要足够长，以便有足够的时间，达到稳定的非平衡态，即热传递达到恒定
• 步长：温度很低，因此步长可以设置的较长
• Initial Temperature：设置为热源、热漏的均值，热源、热漏后面分别设置为 90K、70K，因此此处为 80K

分别设置热源、热漏的温度系综：

注意上图中Atoms in Region，选择了对应的Region。

保存并运行作业。

SCM → Movie 选中默认的能量曲线 del 键删除 → MD Properties → NHCTstat2Energy：

双击坐标轴，选择右侧 Analysis 栏，再选择 LineaRegression 栏，对该曲线进行线性回归拟合： 由于笔者 Movie → Graph → Try To Use Time On E-axes 勾选了，因此曲线的横坐标一般会是时间而非帧数。由于前面的很长时间开始体系热传递还没有达到平衡，因此应该等到足够的时间之后，曲线斜率稳定之后，拟合斜率稳定的这部分数据，因此我们这里作为示意，选择从 10 万 fs 后的数据开始拟合（因为总的模拟时间是 20 万 fs，因此拟合范围为 10 万 fs），然后点击 OK，即得到拟合直线与原始曲线合并显示：

从上图中可以看到，拟合的时候，忽略了前 10 万 fs 的数据。直线数据的标题（上图右上角），可以看到直线斜率为 5.39×10${-5}$ eV/fs(单位为纵坐标单位 ÷ 横坐标单位)

如上文所说，取其一半即 dE/dt = 2.695 ×10${-5}$eV/fs

Model → Lattice，下图中可见，对本例而言，S = a × c = 281.96919 Å$^2$

选中对应的两个原子，即可显示两个原子之间的距离。

带入文中开头的公式中，换算为国际单位制即可求得热导率。

• 模型太小，则热涨落很大，模型越大，结果越趋于稳定。
• 力场的准确度对结果影响很大。
• 该方法对力场、DFTB、MOPAC、DFT均适用。