计算化学公社

标题: 利用Lammps ReaxFF研究反应动力学一例 [打印本页]

作者
Author: Graphite 时间: 2023-7-21 13:20
标题: 利用Lammps ReaxFF研究反应动力学一例
本帖最后由 Graphite 于 2025-4-16 07:48 编辑

利用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时的结构，可以发现确实发生了热解，分子得到了一定程度的聚合，并且得到了苯、长链不饱和烃等物质。（下图左）
当体系中含有氧，即初始投入的物质中包括一些含氧小分子时，可以显著加速这个过程。（下图右）
(, 下载次数 Times of downloads: 167) (, 下载次数 Times of downloads: 161)
RDF也能够以更为定量的方式佐证这一点，如下图所示。
(, 下载次数 Times of downloads: 160)
将含氧体系的pyrolysis.species.out扔进脚本reax_species.py（http://bbs.keinsci.com/thread-38781-1-1.html）中，清洗成表格，作图。此处是对该脚本进行了一定修改，统计的是某类物质的总碳原子数而不是分子数。
可见，在这一过程中，C4-C7、C8-C11的碳数都是先增加再减少，说明是反应的中间产物。C16逐步增加，而且有两个拐点。氢自由基总是保持较稳定水平，说明在链式反应中，氢自由基在生成和消耗中存在动态平衡。这些都指示背后有一定的动力学因素控制。
(, 下载次数 Times of downloads: 151)

pyrolysis.bonds.out可以扔进ReacNetGen里分析反应网络（https://reacnetgenerator.njzjz.win/），下图是前200 ps的ReacNetGen所得反应网络图。
(, 下载次数 Times of downloads: 159)
这里的对应关系：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越多。
(, 下载次数 Times of downloads: 162)
（空间中的电荷分布）
图片有一种无序中带着有序的美感，局部的微结构和氢的转移可能是这个反应的重要因素。

3 总结

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

4 相关工具

现在笔者开发的工具reax_tools也能生成反应的图，支持超大体系、速度很快，见http://bbs.keinsci.com/thread-38781-1-1.html
(, 下载次数 Times of downloads: 90)
(, 下载次数 Times of downloads: 91)

作者
Author: MercuryLamp 时间: 2023-7-23 00:04
本帖最后由 MercuryLamp 于 2023-7-23 12:30 编辑

老师您好，我最近碰到了您之前一个贴中（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一样）

作者
Author: Graphite 时间: 2023-7-23 18:08

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内存还是很有优势的。

作者
Author: MercuryLamp 时间: 2023-7-23 22:14
本帖最后由 MercuryLamp 于 2023-7-23 22:16 编辑

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）

作者
Author: Graphite 时间: 2023-7-24 10:17
本帖最后由 Graphite 于 2023-7-24 10:21 编辑

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

您提到

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

作者
Author: MercuryLamp 时间: 2023-7-24 12:56

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

好的，谢谢您的解答

作者
Author: lmch 时间: 2023-7-29 02:36
请问一下，使用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

作者
Author: Graphite 时间: 2023-7-31 17:52

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

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

作者
Author: lmch 时间: 2023-7-31 23:12

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

感谢指教

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

作者
Author: Graphite 时间: 2023-8-28 09:44

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

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

作者
Author: 2126029128 时间: 2023-8-28 10:02

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

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

作者
Author: ddddgc 时间: 2023-10-22 17:13
本帖最后由 ddddgc 于 2023-10-22 19:06 编辑

老师您好，我对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
......

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

作者
Author: Graphite 时间: 2023-11-7 20:26

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

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

作者
Author: qzupc 时间: 2023-11-8 09:22

Graphite 发表于 2023-11-7 20:26
这个没有统一标准或者统一方法，可以自己调。

好的老师，我调整一下试试

作者
Author: lihaha 时间: 2023-12-27 21:21
答主您好，参考lammps官方文档，感觉应该是reax/c/species会覆盖接邻列表更新的默认设定，而不是reax/c/bonds。不知道我的理解对不对

作者
Author: ggshining 时间: 2024-1-16 11:48
感谢分享，难得的reaxff经验。
“ReaxFF模块可以用KOKKOS的方式使用GPU加速，前提是自行编译了GPU、KOKKOS包”
想请问下，GPU和KOKKOS包都编译了的话，是用哪个包进行加速呢？一起用？咋提交呢

作者
Author: Graphite 时间: 2024-1-17 08:34

ggshining 发表于 2024-1-16 11:48
感谢分享，难得的reaxff经验。
“ReaxFF模块可以用KOKKOS的方式使用GPU加速，前提是自行编译了GPU、KOKKOS ...

前面没写清楚，编译一个即可，一般用kokkos更快
单卡单进程：lmp -in [in文件] -k on g 1 -sf kk -pk kokkkos newton on neigh half
多卡需结合mpirun，详见lammps官网

作者
Author: evenn 时间: 2024-2-1 11:55
老师您好，请问您在文中提到的查看轨迹文件使用的是什么软件呢（我看您的轨迹文件里面有显示键）？我用可视化软件ovito观察输出的dump文件，原子之间没有键连接，而且看到生成的物种和lammps直接输出的species文件所给出的产物类型不同，这个您清楚是什么原因吗？

作者
Author: Graphite 时间: 2024-2-1 15:09

evenn 发表于 2024-2-1 11:55
老师您好，请问您在文中提到的查看轨迹文件使用的是什么软件呢（我看您的轨迹文件里面有显示键）？我用可视 ...

OVITO插件create bonds，根据vdw半径判断成键；OVITO插件load trajectories也可以直接读取lammps fix reax/c/bonds输出的bond文件。
两个单论效果差不多。

作者
Author: Graphite 时间: 2024-2-1 15:14
本帖最后由 Graphite 于 2024-2-1 15:16 编辑

evenn 发表于 2024-2-1 11:55
老师您好，请问您在文中提到的查看轨迹文件使用的是什么软件呢（我看您的轨迹文件里面有显示键）？我用可视 ...

“生成的物种和lammps直接输出的species文件所给出的产物类型不同“，不知道怎么什么生成的物种。
但本质上，软件的算法无非是按照键级/键长的判据把键连关系弄出来，然后再根据键连关系判断物种，所以就是个判据、判断的问题。
vmd/ovito之类是根据键长，可以用键长半径相关的设置调；lammps fix reax/c/species是根据键级，其实也间接和键长相关，可以用cutoff关键词调。
默认判据和自己调都有一定主观指定的意味，所以要做讲究一点的话，得有些依据。

作者
Author: evenn 时间: 2024-2-1 18:09

Graphite 发表于 2024-2-1 15:14
“生成的物种和lammps直接输出的species文件所给出的产物类型不同“，不知道怎么什么生成的物种。
但本 ...

好的，我再琢磨琢磨，谢谢老师回答

作者
Author: 被科研狠狠摩擦 时间: 2024-2-29 21:54
老师，你好，我想请问一下chemtrayzer怎么运行呢

作者
Author: davi 时间: 2024-3-1 11:29
老师，您好，我用lammps但不用反应力场，我看了您的文章还是很有收获，有个问题想要请教一下，再跑lammps之前不需要对分子进行结构优化嘛？

作者
Author: Graphite 时间: 2024-3-1 13:22

davi 发表于 2024-3-1 11:29
老师，您好，我用lammps但不用反应力场，我看了您的文章还是很有收获，有个问题想要请教一下，再跑lammps ...

不管什么软件，MD模拟前结构都要基本合理、和力场匹配。例如原子间距离（角度）等不能和力场定义的平衡长度（角度）相差太远，周期性的体系边缘要处理好等。

能量最小化（结构优化）的意义也就是消除结构的不合理性，避免MD的初始几步因为过高的能量/力/速度而崩掉。如果建模的时候质量已经比较好，能合理地跑起来就不是必须能量最小化（结构优化）。

当然还是看具体问题，例如某炸药正常密度是1.8，想做平衡模拟，却输入了个密度2.1的模型，程序上倒是能跑，但初始瞬时力偏大就直接炸掉了，成了压缩起爆模拟。

作者
Author: davi 时间: 2024-3-1 14:24

Graphite 发表于 2024-3-1 13:22
不管什么软件，MD模拟前结构都要基本合理、和力场匹配。例如原子间距离（角度）等不能和力场定义的平衡长 ...

老师，那是否需要比如采用Gaussian对分子进行opt呢，还是直接用Gsview建好模型将对应力场的长度、角度赋给模型呢，感谢老师

作者
Author: Graphite 时间: 2024-3-1 18:00

davi 发表于 2024-3-1 14:24
老师，那是否需要比如采用Gaussian对分子进行opt呢，还是直接用Gsview建好模型将对应力场的长度、角度赋 ...

不用，guassian优化自然是基于量子化学的方法。既然是做分子模拟用lammps本身、模拟用什么势函数就用什么，minimize一下就可以，如果有必要，跑短步长一小会儿的MD使其平衡即可。

作者
Author: davi 时间: 2024-3-1 18:32

Graphite 发表于 2024-3-1 18:00
不用，guassian优化自然是基于量子化学的方法。既然是做分子模拟用lammps本身、模拟用什么势函数就用什么 ...

感谢老师!

作者
Author: maoxinxina 时间: 2024-4-11 10:29
运行报错ERROR on proc 0: The ReaxFF parameter file ../params/ffield.reax.chon.2019.params cannot be opened: No such file or directory (src/REAXFF/reaxff_ffield.cpp:72)

作者
Author: maoxinxina 时间: 2024-4-11 14:30
把文献中的力场文件输入后，又得到了新的报错ERROR on proc 0: step 2: ran out of space on angle_list: top=16000, max=14179 (src/REAXFF/reaxff_valence_angles.cpp:381)

作者
Author: neoje 时间: 2024-4-14 00:45

maoxinxina 发表于 2024-4-11 14:30
把文献中的力场文件输入后，又得到了新的报错ERROR on proc 0: step 2: ran out of space on angle_list: t ...

您好，能请问下这篇文献是哪篇文献吗

作者
Author: maoxinxina 时间: 2024-4-15 08:27

neoje 发表于 2024-4-14 00:45
您好，能请问下这篇文献是哪篇文献吗

在Amber官网里面的力场文件中能找到这篇文献的来源。

作者
Author: neoje 时间: 2024-4-15 12:41

maoxinxina 发表于 2024-4-15 08:27
在Amber官网里面的力场文件中能找到这篇文献的来源。

好的好的，谢谢！

作者
Author: 王五 时间: 2024-6-9 23:58
本帖最后由王五于 2024-6-10 09:30 编辑

老师您好，感谢您的分享，我在重复您模拟的过程中出现了一些问题，不知如何解决，不知您可否不吝赐教。模拟过程中的in文件和结构文件如上传文件所示，问题：1）在模拟的结果中，并未发现C16，甚至C10以上的物质都很少，不知道是哪一步出现了错误？2）在利用reaxff模拟热解乙炔之前是否要先进行动力学平衡？

作者
Author: Graphite 时间: 2024-6-10 09:51

王五发表于 2024-6-9 23:58
老师您好，感谢您的分享，我在重复您模拟的过程中出现了一些问题，不知如何解决，不知您可否不吝赐教。模拟 ...

（1）这个倒没什么错不错的，用的模型、力场文件和我本就不一样。不过data文件里面边界是341A，物质实际却填充在大概<30A的范围内，是有意为之吗？
（2) 和其他模拟的判断方式一样。

作者
Author: 张德胜 时间: 2024-6-10 11:33

Graphite 发表于 2024-6-10 09:51
（1）这个倒没什么错不错的，用的模型、力场文件和我本就不一样。不过data文件里面边界是341A，物质实际 ...

您前面的例子盒子尺寸不就是341.5A么，应该是按照您的数据进行重复的。

作者
Author: Graphite 时间: 2024-6-10 16:49

张德胜发表于 2024-6-10 11:33
您前面的例子盒子尺寸不就是341.5A么，应该是按照您的数据进行重复的。

我那体系是均匀分布的，是个稀疏气相。他这个一开始所有分子就压在不到盒子边界1%的体积里。

作者
Author: Q1ngKl 时间: 2024-7-9 15:07
请问文中提到的关于bonds和species输出可以在跑完拿到轨迹后重跑获取，具体是怎么操作呢？
是指写一个只dump路径信息的in文件跑完，然后再用个输出bonds和species文件的in文件再跑一遍吗？

作者
Author: Graphite 时间: 2024-7-9 22:20
本帖最后由 Graphite 于 2024-7-9 22:24 编辑

Q1ngKl 发表于 2024-7-9 15:07
请问文中提到的关于bonds和species输出可以在跑完拿到轨迹后重跑获取，具体是怎么操作呢？
是指写一个只du ...

对，最后不是run（步数），而是rerun (相应指令）（相应轨迹）
注意轨迹得够密，如500帧或1000帧，然后reax/c/bonds和reax/c/species的频率和轨迹的频率一致。
如果对bonds和species的采样方式有一些要求（比如reax/c/species每100步采样1次，5次取平均，500步输出一次这种），那么最好还是正常跑的时候输出。

如果你的体系不是特别复杂，只是需要得到物种变化规律，可以用我另一个帖子提供的工具计算物种（0.9版本）。目前只是用van der Waals半径决定分子连接性，所以不太适合特别复杂的结构。这个夏天之后会更新到基于图论的分子识别。也可以用reactnetgenerator，它也可以直接读dump轨迹。

作者
Author: Q1ngKl 时间: 2024-7-10 13:10

Graphite 发表于 2024-7-9 22:20
对，最后不是run（步数），而是rerun (相应指令）（相应轨迹）
注意轨迹得够密，如500帧或1000帧，然后r ...

好的，非常感谢您的解答

作者
Author: ddddnight 时间: 2024-9-17 22:52
# 判断类型1-1、1-2、2-2（也就是C-C、C-H、H-H）成键的键级判据设为0.55、0.4、0.55
老师您好，你这个地方写的是C-C成键的键级判断设为了0.55，我在有些文章里面看PP、PE塑料的热解过程中，对于C-C键的键级判断是0.3，这杨合理吗？这个热解键级判据应该怎么设置呢？
(我参考论文中是这么写的：fix 4 all reaxff/species 1 100 1000 species.out element Cu H C C cutoff 3 3 0.3）

作者
Author: Graphite 时间: 2024-9-18 19:38

ddddnight 发表于 2024-9-17 22:52
# 判断类型1-1、1-2、2-2（也就是C-C、C-H、H-H）成键的键级判据设为0.55、0.4、0.55
老师您好，你这个地 ...

本质上是对键的判定严格不严格。很多文献和文章、ppt常用的是0.3，这意味着成键判定很松。但其实据我挖掘的结果，单纯是很早一篇文献用的0.3，然后lammps默认值是0.3，被沿用下来而已。
像gview/vmd/ovito之类的可视化软件，读取xyz/gro这类没有键定义的坐标文件时，会给<1.6 A/<1.7 A的两个碳之间加个键，和这其实是一码事。
这个值可以调，我这里写成0.55只是因为认为我这个问题（加热小分子合成大分子）应该判定严格一些，免得强静电吸附、半成键不成键的状态被当做成键。关键是道理能说得通、对比的对象保持一致。

作者
Author: ddddnight 时间: 2024-9-18 21:02

Graphite 发表于 2024-9-18 19:38
本质上是对键的判定严格不严格。很多文献和文章、ppt常用的是0.3，这意味着成键判定很松。但其实据我挖掘 ...

好的老师，明白了，感谢您的回答！

作者
Author: 筛石灰 时间: 2024-9-20 21:45
请问老师，lammps reaxff力场可以在windows系统下跑吗

作者
Author: Graphite 时间: 2024-9-21 03:52

筛石灰发表于 2024-9-20 21:45
请问老师，lammps reaxff力场可以在windows系统下跑吗

LAMMPS的所有功能都能在windows下编译使用，但是何必呢

作者
Author: QQi 时间: 2024-9-26 14:50
老师您好，请问这个乙炔热解的例子有相关的文献吗，想学习一下

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-10-18 17:10

Graphite 发表于 2023-7-21 13:20:24
# 设置系综，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

你好，我有个基础问题，想确认一下Lammps用Bussi的CSVR热浴跑NVT的分子动力学要怎么设置。

阅读Lammps手册知道fix nvt默认用Nose-Hoover热浴，同时实现控温和时间积分；而fix temp/csvr用CSVR热浴，只有控温功能而没有时间积分，需要按Thermostats页所说联用fix nve做时间积分来让原子动起来。fix nve默认采用velocity Verlet算法，在有PBC且没有修改受力或速度的fix时产生符合微正则系综的轨迹。

按我的理解是应该如下书写指令：

fix 1 all nve
fix 1 all temp/csvr 360 360 100.0 123456789

主要不确定的在于顺序是不是不能颠倒（会变成NVE），以及这个all作为group ID是否需要提前用set group之类的设置。

作者
Author: Graphite 时间: 2024-10-18 17:53

Uus/pMeC6H4-/キ发表于 2024-10-18 17:10
你好，我有个基础问题，想确认一下Lammps用Bussi的CSVR热浴跑NVT的分子动力学要怎么设置。

阅读Lammps ...

你说的没错。程序上fix nve只是时间积分。nvt就是nve又挂了个Nose-Hoover恒温。其实只是lammps开发者钦定了这种简便写法。

你也可以像fix nve + fix temp/berendsen + fix press/berendsen，得到的就是NpT系综。

能不能颠倒我没试过，正常是先写nve。all是保留关键字，在任何时候下都表示所有原子。

作者
Author: gavin1113 时间: 2024-10-24 20:20
您好，我是新手，请问可以发一下“ffield.reax.chon.2019.params”这个力场文件吗

作者
Author: guoqizhi 时间: 2025-1-11 23:08
老师，请问这个in文件为什么没有进行弛豫？

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

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

作者
Author: Graphite 时间: 2025-3-14 09:37

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

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

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

作者
Author: Graphite 时间: 2025-4-12 09:03
本帖最后由 Graphite 于 2025-4-12 09:36 编辑

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

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

作者
Author: TinnKuka 时间: 2025-5-2 17:49
本帖最后由 TinnKuka 于 2025-5-2 18:13 编辑

老师好，请问为什么我跑完后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 就好了。

作者
Author: yjb 时间: 2025-5-9 21:51
运用老师的输入文件，报错
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不知道为什么还是报错。

作者
Author: Graphite 时间: 2025-5-12 12:50

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恒温器就可以。

作者
Author: yjb 时间: 2025-5-12 21:21

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

谢谢老师

作者
Author: love_yy 时间: 2025-5-22 21:41
各位老师好，能否请教以下对于参照上述帖子已经获得的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” ；
谢谢！

作者
Author: Graphite 时间: 2025-5-23 04:21
本帖最后由 Graphite 于 2025-5-23 04:22 编辑

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

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

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

作者
Author: love_yy 时间: 2025-5-23 12:36

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

谢谢老师回复

作者
Author: jiansi 时间: 前天 13:49
老师您好，请问reax_species.py是否有下载渠道呢

欢迎光临计算化学公社 (http://bbs.keinsci.com/)