利用Lammps ReaxFF研究反应动力学一例

Graphite

利用Lammps ReaxFF研究反应动力学一例

石墨 2023-07-21 首发于计算化学公社

交流QQ群561184358

1、前言

ReaxFF反应力场适用于各类化学反应体系的模拟，在成本-精度的天梯上介于经典动力学（> 数万原子）和半经验AIMD（< 数千原子）之间。
ReaxFF最早是van Duin在2000年前后提出，用来研究一些CH/CHO化合物，后来扩展到N、S、Si甚至金属、稀有气体等体系。

该力场在分类上属于多体势，也就是说，不依赖固定的键连关系，而是只要在一定范围内，就计算相互作用。然后更新距离、速度，如果位置合适，就被认为成键和形成化合物。
这一点和量化/第一性很像，其实ReaxFF的力场参数就是从量化拟合出来的。具体形式不在这篇文章讨论范围之内。
ReaxFF运行时初始结构只需要xyz坐标，带上QEQ之后连电荷都不用自己指定。

因此，ReaxFF最适用于原理简单、动态复杂的反应。例如：小分子热解、燃烧、爆炸、表面融合、吸附、简单基团转移等。
这种反应，要用到点化学性质，但又没必要上到半经验甚至量化。反而是体系可能比较大或考察的情况很多，上到更高级别算不动。
ReaxFF正是因为解决这个痛点而广为流传的。

2、Lammps中ReaxFF的使用：以热解乙炔为例

ReaxFF最早是用van Duin本人开发的同名程序（Fortran编写）计算，后来代码被嵌入到一些程序里，包括ADF、Lammps等。
原程序和ADF都不好用。前者是研究性质的，必须称赞它的科学意义，但毕竟离易用通用还有些距离。ADF是收费软件，把很多功能整合了带个图形界面，卖的很贵，而且封闭、功能固定。
目前最好的ReaxFF软件就是Lammps，更新到C语言的reax/c模块之后，效率提升了很多。Lammps的功能极为强大，虽然速度相比于GROMACS并不快，但是各种技巧、操作是所有分子模拟软件里最全、最实用、最灵活的。
现如今80%以上的反应力场工作的文献，都使用Lammps完成，大部分的ReaxFF脚本、分析软件也都能与Lammps联动。
Lammps中使用ReaxFF时，可以动态产出键的信息、物种分布信息，进而导入ChemTrayzer、ReacNetGen等软件中，进行反应网络分析。
本文仅作抛砖引玉，详细信息见：https://docs.lammps.org/Manual.html，Lammps手册极为详尽，甚至接近教科书的级别。

2.1 程序准备

Lammps中使用ReaxFF模块，需要在编译时包含REAXFF包，同时由于要计算动态电荷，还需要QEQ包。
对于新手，有两种方式准备：
1、自己编译，用cmake编译非常方便，参考Lammps英文手册和《WSL2下Kokkos版加速的Lammps的cmake编译》http://bbs.keinsci.com/forum.php?mod=viewthread&tid=36559
2、现成尝鲜，OpenSuse的软件源里有已经编译好的Lammps（2020年版本、含绝大多数包），支持CPU并行。懒人的办法就是直接用Win10下个OpenSuse版的Linux子系统（WSL），然后sudo zypper install lammps。

ReaxFF模块可以用KOKKOS的方式使用GPU加速，前提是自行编译了GPU、KOKKOS包。
GPU加速的效果和体系的大小有关。但是GPU加速有一个好处是，并行效率的烦恼比较小。这个体系模拟后期可能形成团簇，这样CPU域分解出来不平衡，拖累计算速度（不到刚开始的一半）。但是GPU特别是单卡，基本没有这个顾虑。

因为是多体势，运行的速度直接和密度相关，分子越大、越密集，产生的pair就越多，越慢。
使用16核，5000原子左右的体系，0.2 fs的步长，极低密度的气相热解问题能跑到> 20 ns/day。而5000原子高密度凝聚态可能不到1 ns/day。
虽然大多数用的上ReaxFF的问题几ns、几十ns肯定够了，但是反应力场是必然需要钻研调试的，所以需要特别注意模拟设计的合理性，免得白跑。
因此，在性能调优时，需要综合权衡核数、neighbor、balance方式，避免因为原子跑出CPU核心的”包干区”而导致崩溃。


我在autodl上有一个装好了lmp含reaxff相关包的镜像，有需要可以联系我，将你的账号ID（如4a78977897-xxxx-xxx...）发给我，看到后我会将镜像分享。用此镜像拉起V100实例，然后用lmp-kk -in <in文件> 计算即可。（其实就是个快捷方式，后面可继续加其他需要的参数如-var T 3000等，但不用再加kokkos参数了）

2.2 文件准备

ReaxFF力场使用时最多需要四个文件：in、data、ffield、control
in文件：主要输入，和lj、harmonic之类的势函数不同，ReaxFF的参数不写在in文件里，而是另起一个文件，在in文件里调用它；
data文件：初始结构，这里因为是多体势，只需要坐标和电荷，其中电荷又因为QEq会算，留个空位都写0.0就行。对应的atom_style可以是charge或full。
ffield文件：对于力场参数，如果前人已有，建议先用类似的体系的文献中的，自己拟合会很费事。大部分ReaxFF研究如果改动了力场都会在SI里附力场文件，一般都通用。
只是需要注意SI往往是PDF格式，需要手工复制下来调整，并且与Lammps自带的几个ffield.reax.*文件的格式进行对比。确保格式一致，否则就会报错。特别是数字之间有没有粘连、该有的注释行是否数量一致。
control文件：一些控制计算细节的设置。【这会影响计算结果和稳定性，特别是关注成键、形成物质的时候。】控制文件的内容包括用于计算键的接邻原子的cutoff、氢键cutoff，一般来说，用默认值就足够好，如果确需修改，请参考原文谨慎调整，并且确保修改是有意义的；还有一些关于输出的关键词，实测一般没有改的必要。如果都使用默认值，那么在in文件定义reax控制文件的那个位置写NULL。

下面我们以热解乙炔为例，讲讲data文件和in文件怎么准备。

2.3 data文件
这个体系只有单一的低密度气相。因此，初始建模时，直接根据密度，换算成盒子边长和原子数，用packmol建立均匀放置乙炔分子的盒子。
（in.inp）    filetype xyz
    tolerance 4.0
    seed -1

    output in.xyz

    structure ethyne.xyz
        number 1000
        inside box 0. 0. 0. 341.5 341.5 341.5
    end structure

（linux命令行）
packmol < in.inp

得到in.xyz，然后用VMD将其转成Lammps适用的data文件。用atomsk等软件也可以。
（linux命令行）
vmd in.xyz

（进入VMD 控制台）
topo clearbonds
topo writelammpsdata in.data full

这样我们就得到了一个atom_style为full的结构文件。看一眼：
(in.data)
LAMMPS data file. CGCMM style. atom_style full generated by VMD/TopoTools v1.7 on Thu Jul 20 04:04:11 +0800 2023
4000 atoms
0 bonds
0 angles
0 dihedrals
0 impropers
2 atom types
0 bond types
0 angle types
0 dihedral types
0 improper types
0.0 341.5  xlo xhi   （这里是手动修改了边界，因为.xyz文件不保留边界条件，VMD也不会识别）
0.0 341.5  ylo yhi    (也可以在VMD中先输入pbc set {x y z}，之后再输出的结构文件就包含边界条件）
0.0 341.5  zlo zhi

# Pair Coeffs
#
# 1  C
# 2  H

Masses

1 12.010700 # C
2 1.007940 # H

Atoms # full

1 1 1 0.000000 287.039795 211.846375 315.381836 # C
2 1 1 0.000000 286.540985 212.631454 316.148438 # C
3 1 2 0.000000 286.099945 213.325607 316.826233 # H
4 1 2 0.000000 287.480835 211.152222 314.704041 # H
（以下略）

2.4 in文件
（pyrolysis.lmp）
# 笔者的习惯，开头先定义basename，之后用诸如${basename}.log使得输出的文件名都保持一致。

variable inname string "in"
variable basename string "pyrolysis"

# 单位制

units            real

# 读取结构

atom_style        full
read_data        ${inname}.data

# 调用reax/c形式的势函数来计算原子对间受力、能量，并使用默认控制参数（控制文件名处写NULL）
#  这里和之后的reax/c/species因为都是用的最新最快的C代码，所以统一都是reax/c/*，官方手册上仍然写着reaxff/*，但用法接口都一样。
# 调用力场参数，对所有原子对（* *，两个通配符)，都调用如下参数文件，其中，原子类型1对应参数文件里的C，原子类型2对应H。

pair_style            reax/c NULL
pair_coeff            * * ../params/ffield.reax.chon.2019.params C H

# reax/c自带14个可以计算的量供有需要的人输出，这里只摘出来2个举例，分别是键能、静电作用能。具体c_reax[序号]对应什么物理量见官方手册。
compute             reax all pair reax/c
variable eb         equal c_reax[1]
variable eqeq       equal c_reax[14]

# 应用fix qeq/reax计算动态电荷，每50步更新一次，参考原子的距离上下限为0.0分别8.0，精度要求1e-4，参数用reaxff默认。
#  因为qeq涉及到动态原子列表，所以将更新频率设的太小，会有很大的性能开销。同样的距离上下限也没必要设置的太大。

fix reax_qeq            all qeq/reax 50 0.0 8.0 1e-4 reaxff

# 应用fix reax/c/species输出动态物种数量，每100步采样一次，每50次取平均，总共5000次写入一次文件。
#  判断类型1-1、1-2、2-2（也就是C-C、C-H、H-H）成键的键级判据设为0.55、0.4、0.55
#  注意，这里的采样方式和判据需要根据自己体系的情况。不合理时可能会造成指认出来的物种太过碎片化、太过团簇化、反复震荡等。
#  每5000步输出一次物种坐标（不是特别必要，dump出所有坐标就行了）
#  应用fix reax/c/bonds每5000步输出一次键连关系，并压缩。这可以用来给ChemTrayzer分析。这些输出频率都根据自己的需求来。

fix reax_out_species    all reax/c/species 100 50 5000 ${basename}.species.out cutoff 1 1 0.55 cutoff 1 2 0.4 cutoff 2 2 0.55 position 5000 ${basename}.pos
fix reax_out_bonds      all reax/c/bonds 5000 ${basename}.bonds.out.gz

#  接邻列表设定，注意，fix reax/c/bonds会根据其更新频率，覆盖接邻列表更新的默认设定，所以neigh相关的设定得在那之后。
#  对于这个体系，reax/c/bonds每5000步输出一次够了，但是接邻列表也5000次更新一次的话，大概率是有原子要穿出CPU域分解的空间，导致模拟丢失原子崩掉的。
#  相较于势函数的最大作用半径，再延长2.5个距离单位进行搜索接邻列表，空间形状为盒子型。500步更新一次，不延迟、不校对。
#  需要注意，species和bonds的这些统计输出是依赖CPU和IO性能的，当使用典型的GPU-resident模式跑时，过于频繁地输出会造成极大的性能浪费（超过50%）。此时，可以考虑先不输出，全速跑完拿到轨迹后再rerun统计、输出。

neighbor                    2.5 bin
neigh_modify             every 500 delay 0 check no


# 一般输出设定
# 5000步一次屏显，屏显形式自定义，包括步数、温度、总能；之前定义的reax键能、静电作用能；计算速度、剩余时间。
# 50000步保存一次轨迹。
#  保存日志文件。
thermo          5000
thermo_style    custom step temp etotal v_eb v_eqeq spcpu cpuremain
dump            traj all atom 50000 ${basename}.lammpstrj
log             ${basename}.log

# 设置系综，NVE系综仅作为时间积分使用，用fix temp/rescale对温度实施硬控。（不等于CSVR/V-rescale，不一定适合别的体系)
# 根据温度（1373K）设置初始速度和随机数种子。
fix ensemble         all nve
fix vrescale         all temp/rescale 10 1373 1373 1.1 1.0
velocity        all create 1373 114514

# 能量最小化，设置力和能量的收敛限，最大迭代次数
#  时间步长0.2 fs，运行2500万步，即5 ns。
minimize        1e-4 1e-6 1000 1000
timestep        0.2
run                25000000

# 结束时输出结构
write_data      ${basename}.data

2.5 运行后分析
用mpirun -np [核数] lmp -in pyrolysis.lmp运行模拟
（对于GPU运算，采用KOKKOS方案，如对于单卡：lmp -in pyrolysis.lmp -k on g 1 -sf kk -pk kokkkos newton on neigh half）

运行时，输出大量屏显信息；运行结束，得到了一堆文件，包括：
pyrolysis.data                      最终结构
pyrolysis.log                        运行日志
pyrolysis.lammpstrj              轨迹
pyrolysis.species.out             物种随时间变化，可用来扔给ChemTrayzer和ReacNetGen分析反应网络
pyrolysis.bonds.out.gz                键连关系随时间变化，同上
pyrolysis.pos                        物种中心坐标随时间变化，可以自行作图

我们先查看轨迹，观察0 ns和 5 ns时的结构，可以发现确实发生了热解，分子得到了一定程度的聚合，并且得到了苯、长链不饱和烃等物质。（下图左）
当体系中含有氧，即初始投入的物质中包括一些含氧小分子时，可以显著加速这个过程。（下图右）

RDF也能够以更为定量的方式佐证这一点，如下图所示。

将含氧体系的pyrolysis.species.out扔进脚本reax_species.py（http://bbs.keinsci.com/thread-38781-1-1.html）中，清洗成表格，作图。此处是对该脚本进行了一定修改，统计的是某类物质的总碳原子数而不是分子数。
可见，在这一过程中，C4-C7、C8-C11的碳数都是先增加再减少，说明是反应的中间产物。C16逐步增加，而且有两个拐点。氢自由基总是保持较稳定水平，说明在链式反应中，氢自由基在生成和消耗中存在动态平衡。这些都指示背后有一定的动力学因素控制。


pyrolysis.bonds.out可以扔进ReacNetGen里分析反应网络（https://reacnetgenerator.njzjz.win/），下图是前200 ps的ReacNetGen所得反应网络图。

这里的对应关系：1：C；2：CH；3：H2；4：CH4；5：H；6：CH2；7：CH3；8：C2H2；9：C2H。
这说明前200 ps中。乙炔先裂解，然后是大量的C1、C2物种的相互转化，转化过程中，氢占据重要地位。这一点与前面一致。


将pyrolysis.pos进行处理和作图，考虑到物种太多，而H总是偏正电，C总是偏负电，可以用分子电荷的方式体现空间中的物种分布。电荷越正，H越多。

（空间中的电荷分布）
图片有一种无序中带着有序的美感，局部的微结构和氢的转移可能是这个反应的重要因素。

3 总结

ReaxFF反应力场适用于各类化学反应体系的模拟，特别是原理简单、动态复杂的反应/过程。
利用合适的工具和工作方法，我们可以实现一定程度上的“虚拟实验”，得到一些实验室里无法得到的信息，例如燃烧、爆炸、分解反应的历程、连续复杂过程中物质的相互转化等。

4 相关工具

现在笔者开发的工具reax_tools也能生成反应的图，支持超大体系、速度很快，见http://bbs.keinsci.com/thread-38781-1-1.html

jiansi

老师您好，请问reax_species.py是否有下载渠道呢

love_yy

Graphite 发表于 2025-5-23 04:21
直接用reax_tools或者reacnetgenerator读轨迹分析得了，说实话基于LAMMPS键级和cutoff，有时候反而结果相 ...

谢谢老师回复

Graphite

love_yy 发表于 2025-5-22 21:41
各位老师好，能否请教以下对于参照上述帖子已经获得的lammps.trj轨迹如果rerun调整bond_order cutoffs来重 ...

直接用reax_tools或者reacnetgenerator读轨迹分析得了，说实话基于LAMMPS键级和cutoff，有时候反而结果相当的紊乱

reax_tools用vdw缩放因子调意外地又LOW又好用

love_yy

各位老师好，能否请教以下对于参照上述帖子已经获得的lammps.trj轨迹如果rerun调整bond_order cutoffs来重新得到species，在rerun处应该怎么写呢，我根据官网和上述模板在去掉能量最小化和初始速度后对run部分作出了如下修改，一直报错“ERROR: Read_dump must use at least either 'id' or 'type' field (../read_dump.cpp:1175)“，请问正确的方法是什么？
（1）直接将run部分修改为“rerun ${basename}.lammpstrj dump every 1 id type xs ys zs”
（2）将trj文件转化为xyz，将run部分修改为：“rerun ${basename}.xyz dump every 1  x y z” ；
谢谢！

yjb

Graphite 发表于 2025-5-12 12:50
1.可能是文件保存或者同步问题，这里识别成0.0了
2.看fix temp/rescale的手册，window=10.0其实没意义， ...

谢谢 老师

Graphite

yjb 发表于 2025-5-9 21:51
运用老师的输入文件，报错
ERROR: Computed temperature for fix temp/rescale cannot be 0.0 (src/fix_te ...

1.可能是文件保存或者同步问题，这里识别成0.0了
2.看fix temp/rescale的手册，window=10.0其实没意义，而且前面也说了temp/rescale不一定适合所有体系，体系不特殊的话正常的nose-hoover或者csvr恒温器就可以。

yjb

运用老师的输入文件，报错
ERROR: Computed temperature for fix temp/rescale cannot be 0.0 (src/fix_temp_rescale.cpp:136)
输入文件如下：
fix ensemble         all nve
fix vrescale        all temp/rescale 10 600.0 600.0 10.0 1.0
velocity        all create 600.0 114514
明明设置了温度为600不知道为什么还是报错。

TinnKuka

老师好，请问为什么我跑完后rerun dump file的时候fix reaxff/species 输出的是单个原子呢？像这样：
#  Timestep    No_Moles    No_Specs           O           C           O           C           H           C           O           H
       1000        5550           8         500         500         500        1500        2020         510          10          10

我试过更改cutoff结果不变。

------

已解决！我的fix reaxff/species 刚开始Nevery Nrepeat Nfreq设置的是 10 10 1000， 改成 1 1 1000 就好了。

Graphite

Gonglinquan 发表于 2025-4-11 22:31
您好，我最近在使用您开发的reaxff_species.py代码进行物种数据的清洗，这个代码在一个较小的species.out文 ...

0.3版本已经不再维护，请用最新版本

Gonglinquan

您好，我最近在使用您开发的reaxff_species.py代码进行物种数据的清洗，这个代码在一个较小的species.out文件，1M左右，能够正常得到输出结果，但在另一个较大的文件，有10M左右，却只能将所有原子视为一个分子，无法进行详细识别，请问这可能是什么原因造成的呢？

Graphite

菲比的洪均 发表于 2025-3-13 11:33
老师您好，最近在使用LAMMPS的时候遇到了一个问题，就是在计算一个只包含CHN三种元素的体系在高温下的分解 ...

弛豫阶段别加高温，或者初始建模就不合理。

菲比的洪均

老师您好，最近在使用LAMMPS的时候遇到了一个问题，就是在计算一个只包含CHN三种元素的体系在高温下的分解路径时，在弛豫阶段体系就会发生分解，请问您有什么解决办法来防止体系在弛豫阶段的分解吗

lbl

老师您好，我想知道该怎么rerun呢，我在rerun的时候出现了丢失原子的报错，正常跑则没有这个报错，而且我rerun出来的species文件中没有形成团簇，只有单个原子，我想知道该怎么解决呢？

guoqizhi

老师，请问这个in文件为什么没有进行弛豫？