ReaxFF-热导率方法2：T-NEMD方法结合传热功率统计

热导率计算公式：

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。

保存并运行作业。

结果分析

1，dE/dt的求取

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

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