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

Graphite

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

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

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包。
使用单张RTX 3090时，速度和服务器32核差不多; 使用单核+单张Tesla V100时，某些时候甚至不如32核。
GPU加速的效果和体系的大小有关，体系太小，可能还不如10核。但是GPU加速有一个好处是，并行效率的烦恼比较小。这个体系模拟后期可能形成团簇，这样CPU域分解出来不平衡，拖累计算速度（不到刚开始的一半）。但是GPU特别是单卡，基本没有这个顾虑。

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

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

MercuryLamp

老师您好，我最近碰到了您之前一个贴中（http://bbs.keinsci.com/thread-18771-1-1.html）同样的问题，发现自己编译的kokkos加速版本lammps使用多核运行速度反而远远慢于单核运行（我这边配置是7950x+4090），您这个贴子中也提到

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

，请问使用kokkos加速版本的lammps是不是应该单核运行配合GPU加速才最好呢？有没有什么办法能在GPU加速的同时用多核并行进一步加快计算速率呢？

我看您之前的这个贴子中付老师提到

核数越多越慢是非常正常的现象，现在MD引擎中绝大部分计算都是GPU完成，用多个CPU核心大多数时候的作用就是增加通讯消耗

刘老师也提到

直接mpirun -np  4 好像不太行，你看一下是不是这样GPU也并行了，cpu和gpu上都会分4个任务并行，速度会减慢。
要设置1个gpu任务调用多个CPU线程吧？

不知lammps能否像刘老师所说的只设置一个GPU任务调用多个CPU线程。（类似于gromacs一样）

Graphite

MercuryLamp 发表于 2023-7-23 00:04
老师您好，我最近碰到了您之前一个贴中（http://bbs.keinsci.com/thread-18771-1-1.html）同样的问题，发现 ...

您好。关于几核带一卡，目前我试的结果是，用MPI带多进程，会开出来一摸一样的几个一核带一卡的任务，互相无意义地抢资源。用openMP带多线程，速度与单线程持平。说明Lammps在使用gpu加速时几乎是把所有负载都推到gpu上了。从原理上说不是所有计算都适合用gpu算，比如gromacs就是一部分cpu一部分gpu效果最好。这可能也是一个lammps在很多gpu加速中表现一般的原因。
目前我的结论是lammps确实GPU提升一般，对于单机和中等规模集群性价比都很暧昧，达不到gromacs那种几乎必入gpu的情况。比如万元能买的gpu跟万元能买到的cpu+内存互相都拉不开两三倍的差距。
其实lammps对一些常见的势函数，gpu加速还是可以的，只不过那些体系还不如用gromacs做更快。
我最近也在研究工作站和小集群lammps性能调优这个事，过段时间可能出个系统的测试报告。
对于您的机子，可能现阶段最好的办法是开两个任务，一种一个单核带单卡，另一个纯cpu运行，7950x搭ddr5内存还是很有优势的。

MercuryLamp

Graphite 发表于 2023-7-23 18:08
您好。关于几核带一卡，目前我试的结果是，用MPI带多进程，会开出来一摸一样的几个一核带一卡的任务，互 ...

好的，非常感谢老师您的建议，还有一个问题想请教一下您

您提到

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

请问这个意思是指in文件中写reax/c和reaxff是等效的吗？因为我在看lammps手册的时候似乎并没有找到reax/c关键词，只看到了reaxff，同时google的时候看到似乎在最近的版本中已经将reax/c重命名为reaxff了（参考：https://matsci.org/t/pair-style- ... ammps-29oct20/40964）

Graphite

MercuryLamp 发表于 2023-7-23 22:14
好的，非常感谢老师您的建议，还有一个问题想请教一下您

您提到

我用的版本是23jun22，这个版本上以reax/c为主，对我的体系也更快。后期有没有重新澄清命名我不太清楚，确实是有点乱的，可以都试试。

MercuryLamp

Graphite 发表于 2023-7-24 10:17
我用的版本是23jun22，这个版本上以reax/c为主，对我的体系也更快。后期有没有重新澄清命名我不太清楚， ...

好的，谢谢您的解答

lmch

请问一下，使用kokkos包模拟，运行命令中关键词 newton on neigh half 对计算结果有影响吗？此外，comm no这个关键词对计算结果有影响吗？
如果仅仅是对计算速度/性能有影响，那就不用过于考虑这些关键词的设置了。
附官网解释
neigh value =half  half neighbor list built in thread-safe manner
newton =  on set Newton pairwise and bonded flags on
comm/reverse value = no  perform communication pack/unpack in non-KOKKOS mode

Graphite

lmch 发表于 2023-7-29 02:36
请问一下，使用kokkos包模拟，运行命令中关键词 newton on neigh half 对计算结果有影响吗？此外，comm no ...

原理上都没什么太大影响，只是neigh和newton在某些GPU平台上，不加报错而已。

lmch

Graphite 发表于 2023-7-31 17:52
原理上都没什么太大影响，只是neigh和newton在某些GPU平台上，不加报错而已。

感谢指教

2126029128

老师你好，请问那个反应网络中的数字123，这些内容是怎么定义的呢

Graphite

2126029128 发表于 2023-8-27 18:03
老师你好，请问那个反应网络中的数字123，这些内容是怎么定义的呢

reacnetgen这个软件自己动态定义的，它会附各个序号对应的分子式。不过这个分子式里面体现的单键双键只是估的，没有明确化学意义，需要注意。

2126029128

Graphite 发表于 2023-8-28 09:44
reacnetgen这个软件自己动态定义的，它会附各个序号对应的分子式。不过这个分子式里面体现的单键双键只是 ...

好的，谢谢老师的回复。老师，这个定义是moname文件中的123吗？另外利用.pos文件来进行绘制电荷分布图，是如何实现的呢

ddddgc

老师您好，我对bonds文件有些疑问，输出的结果中单键（nb=1）的原子，mol都为0，这是不是错误的输出？这里的分子id是如何定义的，望老师不吝赐教

这是in文件中主要的关于ReaxFF部分的描述
pair_style      reax/c control.reax lgvdw yes checkqeq yes safezone 1.6 mincap 100
pair_coeff      * * ffield.reax.lg C H O N

fix  bond2   all  reax/c/bonds  1000 bonds.reax        
fix  1  all qeq/reax 1 0.0 10.0 1e-6 reax/c


其中control.reax文件内容如下
simulation_name         enery_material  ! output files will carry this name + their specific extension

tabulate_long_range     0 ! denotes the granularity of long range tabulation, 0 means no tabulation
energy_update_freq      0
remove_CoM_vel          500 ! remove the trans. & rot. vel around the CoM every 'this many' steps

nbrhood_cutoff          4.5  ! near neighbors cutoff for bond calculations in A
hbond_cutoff            6.0  ! cutoff distance for hydrogen bond interactions
bond_graph_cutoff       0.3  ! bond strength cutoff for bond graphs
thb_cutoff              0.001 ! cutoff value for three body interactions

geo_format              0    ! 0: xyz, 1: pdb, 2: bgf
write_freq              5000   ! write trajectory after so many steps
traj_compress           0    ! 0: no compression  1: uses zlib to compress trajectory output
traj_title              dump ! (no white spaces)
atom_info               0    ! 0: no atom info, 1: print basic atom info in the trajectory file
atom_forces             0    ! 0: basic atom format, 1: print force on each atom in the trajectory file
atom_velocities         0    ! 0: basic atom format, 1: print the velocity of each atom in the trajectory file
bond_info               0    ! 0: do not print bonds, 1: print bonds in the trajectory file
angle_info              0    ! 0: do not print angles, 1: print angles in the trajectory file


下面是data文件的部分说明
43776 atoms
4 atom types
0.0 68.0212 xlo xhi
0.0 78.7312997002 ylo yhi
0.0 75.3597982121 zlo zhi

Masses

1 12.0107  # C
2 1.00794  # H
3 15.9994  # O
4 14.0067  # N

Atoms  # charge

1 1 0.0 63.9696728539 39.2140006126 1.6126996817
2 1 0.0 67.5503518095 30.9267440061 1.6126996817
......

qzupc

老师您好，针对您提到的输出动态物种数量要增加各原子间成键的键级判据，即cut off之后的阈值设置，我想请问下这个针对自己的体系的原子键级判据是通过什么可以查或者计算得到呢？
我目前的体系情况是我所加的K2CO3被识别成了很多其他碎片化的物种类型，考虑可能是由于我没有加这个原子成键判据导致的，我之前的in文件在species一行中没有设置cut off阈值

Graphite

qzupc 发表于 2023-11-6 22:55
老师您好，针对您提到的输出动态物种数量要增加各原子间成键的键级判据，即cut off之后的阈值设置，我想请 ...

这个没有统一标准或者统一方法，可以自己调。