REAXFF的性能调优与一点经验 - 计算化学公社

1、前言
自从前帖《利用Lammps ReaxFF研究反应动力学一例》、《LAMMPS中密度不均一体系域分解问题，和性能调优的一点经验》、《Lammps reaxff物种统计分析工具》发布以来，受到大家的广泛关注，不胜荣幸。在模拟实践中，大家也常常遇到许多问题，有些具有共性，在这里愿共同探讨、抛砖引玉。

ReaxFF本身是效率与精度间的妥协、是分子模拟的一个主要矛盾的冲突点和平衡点。分子模拟处处是矛盾。我们无法消除固有矛盾，只能谦卑地调和矛盾。有一天（可能用不了几年）ReaxFF会被机器学习势、半经验AIMD等等取代，进入矛盾的下一个阶段，不过我们现在仍然有讨论ReaxFF的必要。

2、LAMMPS ReaxFF计算的方式和流程
分子动力学模拟无非是原子小球受周围环境影响进行运动。根据运动方程和相互作用函数（力场），对运动进行演绎。这样，模拟大致的流程如下：

（1）确定原子周围可能有影响的其他原子
（2）基于周围环境，为原子计算受力、加速度，更新速度和位置
（3）每过一段时间，计算和输出一些我们感兴趣的量

其中步骤（1）称为邻居搜索，它的范围我们姑且称作”感知范围”，这一步往往很耗时。步骤（2）也就是我们熟知的受力计算和动力学走步，我们希望性能全都在这里释放。步骤（3）也会造成一定耗时。

3、健壮性
模拟的健壮性主要受3部分影响：
（1）结构的合理性、力场适配性、外界附加条件（温度、外力等）的强度。
（2）动力学走步的步长，它决定模拟的精细度，步长越短，运动也就越平缓、模拟越健壮。
（3）“感知范围”，由两部分组成：（1）键级计算范围，也就是原子考虑与它的邻居发生成键、断键、干涉其他键的范围，这个在control文件中，对有机物为主的体系，默认的4.5 A是个还不错的距离，对于一些无机物，可以考虑加大。键级计算范围对模拟的精度和结果影响很大。（2）缓冲范围，也就是在键级计算范围以外预留的缓冲区的大小，即neigh设定的距离。缓冲越保守（缓冲范围大、缓冲重新计算频率快），越不容易因为瞬时受力、瞬时高速使模拟出错。

4、性能
影响性能的因素有4点：
（1）步长，显而易见。
（2）模拟的密度，密度越大，原子感知范围内的其他原子越多，性能开销越大。
（3）感知范围越大，性能开销显著增大。
（4）额外计算和文件输出的开销。

5、健壮性与性能的妥协
步长直接影响模拟精度，需优先确定，对于温和缓慢的反应，步长可以设为0.25-0.5 fs，高温、大外力、不平衡的体系需要低至0.1 fs。密度也是我们无法决定的。因此，性能调优的关键在于感知范围的调整。
感知范围越小，那么原子的运动越“盲目”，自然模拟不健壮。感知范围太大了，凭空增加开销——一般来说6-8 A以外的原子对当前原子的化学性质影响很小，而且也不会毫无征兆地“侵入当前原子的领地”造成异常。

6、调优策略
（1）确定当前体系适用的步长
（2）根据体系的性质（结构本身是否平衡、温度、外力的大小）确定control文件中键级计算范围，和in文件中neigh的半径和频率。这里我们姑且将短半径、低频率称为“激进”（向健壮性作了妥协，换取了性能），反之则为“保守”。
（3）确定额外计算和文件输出的频率。由于邻居搜索和ReaxFF的物种输出(reax/c/species)、键级输出(reax/c/bonds)都是CPU、内存为瓶颈的计算，而且reax/c/species等有时还会覆盖neigh频率，尽量使后者为邻居搜索频率的整数倍，以使得任务周期同步。例如neigh频率为每50步计算一次，reax/c/species可设为100 5 500（100步计算1次，5次取平均，500步输出一次），轨迹输出为1000步输出一次。

7、上中下三策
（1）激进策略：温度不高/本身反应性不强/模拟已经大差不差开始堆时间，步长0.25-0.5，neigh频率50-100，半径1.5-2.0。
（2）温和策略：步长0.2-0.25，neigh频率20-50，半径2.0-2.5。
（3）保守策略：温度很高/反应性强/不平衡性强/初始结构不太好稍微跑会儿放松下，步长0.1，neigh频率10、半径>=3.0。

8、特殊情况
如果保守策略也跑不了，那90%概率是初始结构不好、力场不适配、条件偏极端等，一般来说3000 K或10 km/s MSST冲击使含能材料爆炸，保守策略也是能正确跑的。此时应该：
（1）检查力场的适用性，力场有的为高温/低温开发，有为稳定聚合物/不稳定的含能材料开发等，假如温和的反应用了含能材料的力场，那必然招致意外的分解；
（2）结构的合理性，如果结构还能救，可以先能量最小化，再用非常保守的策略一点点放松，如0.005 fs跑5000步、0.02 fs跑5000步...渐渐增至0.1 fs。
（3）外力和边界条件施加的合理性，例如外界条件先不要立刻施加，平稳后再一点点施加。另外像unit metal下，输入速度的单位是A / ps，力的单位是eV / A，有时会弄错。周期性重复基底的边界处理、边界有反射壁或LJ壁的也要小心。

9、GPU计算时的设置
由于目前版本LAMMPS的限制，截止2024年初，一般性的解决方案是用单张如Tesla V100一类的，具有双精度浮点的显卡进行计算，CPU使用单进程。LAMMPS需要正确安装KOKKOS模块。此时的运行指令为lmp -in [in文件] -k on g 1 -sf kk -pk kokkos newton on neigh half。1.0密度下，激进策略通常能达到100万原子步/秒。
要试用Tesla V100，可在autodl平台上租几个小时。二手显卡也很便宜。