ReaxFF-热导率方法1：T-NEMD方法结合温度分布

计算方法与思路

在体系中划分一个区域为“热源”，类似划分一个区域为“热漏”，二者之间存在一个热传导区域。系统不断通过“热源”以一定速率释放能量，而“热漏”则以相同速率吸收能量，在经过足够时间长度之后，整个系统的温度将达到一种平衡状态。热传导区域的温度分布曲线，将会存在一个稳定的梯度：靠近“热源”的位置温度最高，靠近“热漏”的位置温度最低。

温度分布曲线：

A、B、C三个方向（如何显示方向？），分别均分为N个切片，温度曲线是这N个切片的温度。会得到三个温度分布函数，对应三个方向的温度分布。每个切片的温度，是这个切片内所有原子，从上一帧到当前帧的时间范围内的平均温度。因此切片数不建议太多（以确保片内原子个数较多，具有统计平均的意义），保存轨迹的频率不宜高（默认100建议改为1000或2000甚至根据需要改到更大数值），否则都会导致温度曲线震荡会比较剧烈，失去统计意义。

热导率计算公式：

k = W/(S⋅grad(T))

其中，

W为释放能量的速率，软件中的单位为Hartree/fs，换算为国际单位：1 Hartree/fs = 4.35952 × 10^-3 W
S为热传导方向的横截面，软件中的长度单位为Å，面积为Å²，换算为国际单位：1 Å² = 1 × 10^-20 m²
grad(T)是温度曲线的梯度，软件中得到数据单位是K/Å，换算为国际单位：1 K/Å = 1 × 10¹⁰ K/m

计算引擎：

本文以ReaxFF为例进行演示，实际上使用普通力场的分子动力学、DFT、DFTB分子动力学均可用于热导率计算，除了Main面板（引擎自身参数）的设置不同，其他设置是完全一致的。力场的质量会直接影响热导率计算的准确性。

关于周期性：热导率计算对于0~3维周期边界条件均支持，用户需要注意在没有周期性的情况下，哪个方向，哪一段区域是热传导区域，找到对应的温度分布数据即可。

软件版本要求：

AMS2022.101或以上

模型与参数

AMS生成的温度曲线是从坐标原点开始的，为了和模型的位置完全对应，我们最好将Cell的顶点定义为坐标原点，如下图所示：正常设置分子动力学基本参数。需要注意的是，为了确保温度曲线数据稳定，我们将保存轨迹的频率从默认的100改为了3000，总共运行60万步。关于步数：需要让体系的温度分布达到平衡状态，多少步才能达到平衡？模拟之前，我们是不知道的，因此可以设置的大一些。设置温度：如果我们关心什么温度下的热导率，这里就设置为多少，这个温度将会是体系的总体平均温度。分别设置热源区域（Source）和热漏区域（Sink）（如何添加Region，参考教程如何创建分区）实际上这里只有Souce和Sink用得上，两个传热区也顺便设置（注意传热区2虽然看起来是分成两段的，但由于存在周期性，因此两端实际是连在一起的）。为了确定两个传热区的长度一样，从而热源沿着两个方向传导的功率是一样的（详细参考下文中关于W的设置），点击Region名字后面的√，可以选中该Region所有原子，从而在窗口底部看到原子数目，确认传热区1、传热区2原子个数一样。

设置热传导模型：上述参数：

Method是实现热传导的三种算法：1) Simple方法中，原子速度按比例放大（或缩小）一个因子，该因子与能量速率对应，不考虑总动量的守恒；2) HEX 方法中，速度以总动量守恒的方式缩放，这会在总能量中引入一个小但可见的漂移；3) e HEX 算法，通过添加三阶时间积分校正来消除漂移，从而改进HEX。我们测试SiO₂的热导率，Simple与HEX几乎没有差别，e Hex略有小幅变化（仅1‰量级）。
设置为2000步以后才开始进行热传导（因为MD前面的一些step，体系有可能还没有达到平衡，这个过程相当于让体系充分弛豫，这个值也可以适当调整
能量传递速率（即功率）为0.004 Hatree/fs，对应热导率公式中的W。注意，W是指热源传递出去热量的总功率，本例中的周期性原因，它向±C方向都有传递，因此在其中一个方向上，只有一半的功率。基于这一点而言，两个热传递区域设置为一样长更严谨，能够确保功率为总功率的一半，即公式中W的值为Heat Rate值的一半。该值的设置太大，有可能过载（系统承受不了这么大能量流，计算会报错：Error: Trying to pump Too much heat from the system. Decrease the HeatingRate value），太小则达到平衡的速度太慢，需要的总MD步数会很高。
指定热源为Source区域的原子，热漏为Sink区域的原子。当然用户也可以不使用region，直接勾选Show，然后在下面对应的ABC数据中设定热源与热漏的范围，ABC为晶格矢量，后面的数据为分数坐标，用户可以通过修改数值，观察阴影区域的变化，从而正确理解数据的作用。

保存作业，并修改生成的*.run文件（记事本可以修改），在其中的Trajectory字段增加一行TProfileGridPoints 30：

    Trajectory
        SamplingFreq 3000
        TProfileGridPoints 30
    End

表示体系沿着ABC三个方向，分别得到30个温度数据。

保存*.run文件后，AMSJobs 选中该作业，Job → Run运行。

温度曲线生成

计算完毕后，SCM → Kf Browser → File → Expert Mode，窗口底部输入tempprofile，在MDHistory中可以找到每一帧（本例60万步，每3000步保存一次，加上初始结构因此有201帧）的温度分布函数数据，例如101帧C方向：

这就是温度曲线纵坐标值。

数据处理小技巧:

将所有数值复制到Word文档中，选中数据之间的间隔，用查找替换功能，将其替换为换行符^p，再选中剩余的不规则字符，类似替换为空字符（什么也没有），则变成了规则的一列数据。

横坐标值就是从0~C之间的等分30个数值，如下图所示，C=267.448384 Å，每个数据点的横坐标间隔为C/30，热导率公式需要的C方向的横截面积S为4002.8113 Å²（即热导率公式中的S）：

数据整理后作图（为了确保已经达到热力学平均，这里我们列出了最后一帧即60万步以及30万步的数据，二者基本上重叠，因此可以认为30万步的时候已经达到了热力学平均，否则用户需要设置更长的MD steps）：

图中数据横坐标单位为Å，纵坐标单位为K，温度最高处为热源，温度最低处为热漏，因此温度有上升、下降两段
拟合出温度梯度，单位为K/Å
梯度取绝对值即热导率公式中的grad(T)（因为源和漏的交换位置，则梯度反号），拟合上图曲线的梯度，得到约为3.9 K/Å。

至此热导率公式中所有数据均齐备了，用户可以自行带入公式计算热导率：

k = W/(S⋅grad(T))

= (0.002 Hartree/fs) / (4002.8113 Å² * 3.9 K/Å )

= 8.719 * 10^-6 W / (4.0028113*10^-17 m² * 3.9 * 10¹⁰ K/m)

= 5.59 W/mK

注意：

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

目录

ReaxFF-热导率方法1：T-NEMD方法结合温度分布

计算方法与思路

模型与参数

温度曲线生成

注意：