计算化学公社

标题: BaneTop：Gaussian力场参数转换器 [打印本页]

作者
Author: wal 时间: 2025-5-20 22:28
标题: BaneTop：Gaussian力场参数转换器
本帖最后由 wal 于 2025-5-21 15:08 编辑

5.21 小更新，修了一个罕见情况下filesystem转换不了绝对路径报错的小bug。没遇到的不更也无所谓，没有动主逻辑

发帖的时候犹豫了一会，这个小玩具本身的功能全是处理力场参数的，发QC板块好像有点奇怪，但是在动力学板块发把rtp文件转换成gjf的帖子好像更奇怪（
最后因为本程序的目的是服务于Gaussian做ONIOM(QM:MM)，所以还是发在这里了，如果sob老师觉得不合适求帮忙移过去

如前所述，本程序的目的是解决ONIOM计算中配体爆缺参数的情况。笔者去年9月左右努力过ONIOM(QM:MM)，结果一上来就被塞满一整张屏幕的缺参数报错吓晕。不仅如此，硼这玩意我还没还没查到现成的文献给参数，于是直接放弃了。然而最近做的课题大概是不得不考虑完整蛋白质了，所幸在此期间见识长进了不少，了解到了sobtop程序，于是与Claude3.7、Gemini2.5等商议了一阵，摸了个小玩具出来。
本程序使用C++17编写，数学计算部分由lapack完成，在ubuntu22.04中使用gcc11.4编译。sob老师没有公开sobtop的源代码，故本程序将借助模板文件进行交互。由于ubuntu22.04的glibc库比较新，笔者测试多数旧服务器出现兼容性问题，故采用完全静态链接(除了一个啥玩意系统库只有动态版来着我忘了)，在rocky8.10服务器测试可以正常使用。
BaneTop写成依托的源码暂时不公开了，不然你们就会发现笔者其实全靠AI，自己根本不会C++()

准备工作
I.按照sobtop主页的提示安装sobtop。如果你要算的配体含有不常见原子，你需要在assign_AT.dat里写好自定义辨认规则，并在LJ_param.dat里添加对应原子LJ势参数。例如，笔者算的体系是BODIPY，硼连4个键，于是在assign_AT.dat中添加：

$B_NNFF
nbond 4

复制代码

在LJ_param.dat里添加：

B_NNFF 1.772 0.095 ;Ref:10.1049/mnl.2009.0112

复制代码

II.下载压缩包 (, 下载次数 Times of downloads: 9) 并解压到一处，修改conf文件里的sobtop.path，使之指向sobtop可执行文件，还要修改sobtop.templates_dir，指向压缩包里自带的templates

# Sobtop可执行文件
sobtop.path = /path/to/sobtop
# Sobtop模板目录
sobtop.templates_dir = /path/to/templates

复制代码

完成后，将这个文件移动到~/.banetop，这样程序就可以直接加载了。随后添加可执行权限并进行恰当设置使得banetop可以直接输入命令调用：

chmod +x banetop
echo 'alias banetop=path/to/banetop' >> ~/.bashrc
# 或者加进path
echo 'export PATH=path/to/banetop:$PATH' >> ~/.bashrc

复制代码

文件方面把做了freq的fchk文件、与此fchk原子顺序一致的mol2文件(、包含RESP电荷的chg文件)放在一处即可。看这个红字，你猜我把苯环崩碎多少次才发现顺序不一致会导致力参数瞎配把分子模拟崩

1 不涉及周期势二面角参数缺失
I. Gaussian辨认的优先，认不出来的求助sobtop
由于sobtop和Gaussian在原子类型的指认情况有所不同，参数值可能也因使用的版本原因有所差异(e.g. 对于CA-CA-HA，g16 A.03自带的谐振弯曲力常数是35，sobtop里的是50)
这种情况需要先随便交一次amber计算，得知Gaussian究竟知道多少参数。
在gview里把计算方法改成mechanics，Amber，然后保存提交
(, 下载次数 Times of downloads: 32)
得到如下满满一屏缺参数报错
(, 下载次数 Times of downloads: 33)
(Gaussian你知道的是不是有点太少了！)
得到这个log后，与fchk文件放在一起，运行：

banetop -s genpara filename.mol2 filename.fchk
banetop -u filename.log param.txt

复制代码

运行完毕后，屏幕上应当打印所有缺失的力场参数，这些参数都是sobtop基于mSeminario方法计算的。
(, 下载次数 Times of downloads: 29)

把这些参数空一行复制到gjf文件末尾，写上Amber=softfirst，再次提交计算就能正常结束了。(图略)

如果想要用RESP电荷，把chg文件也放在当前目录下，执行：

banetop -c molecule.gjf charges.chg

复制代码

这样就把RESP电荷写进gjf文件了。

II. 先进的Sobtop已经完全战胜了落后的Gaussian！
由于Sobtop非常强大，我们完全可以跳过Gaussian，直接使用sobtop产生的参数：

banetop -g fchk filename.fchk

复制代码

使用此命令时，程序会依次做这些事：

检查同一目录下有没有同名mol2，如果没有则尝试调用obabel命令行工具进行转换；如果都不可用，报错退出
调用sobtop产生包含全部所需信息的rtp文件，生成参数的策略是指认Amber原子类型，尽量使用原力场参数，缺失的用mSeminario方法计算。如果对此不满意可以自行修改templates目录下的genrtp.txt，这个文件在运行时替换一下占位符就重定向给sobtop了，可以按照按照sobtop操作步骤改成自己喜欢的样子，比如用gaff之类的。
将rtp文件里的所有参数进行单位转换，并翻译成Gaussian能看懂的格式(哇去Gaussian的谐振弯曲/拉伸势函数跟gromacs居然不一样，gromacs的比Gaussian的多带个1/2系数，这个邪门bug困扰了我好几天)
使用sobtop导出带键连关系的gjf，把原子类型、力场参数全都写进去(此处如果同一目录下有同名chg文件会把电荷也一起写进去)

这样导出的gjf文件是可以直接使用amber=(softfirst,firsteq)开始计算的，不会报任何缺参

2 从谐振势拟合周期势二面角参数
写本程序期间，在查Gaussian说明书的时候，发现了一个很头疼的问题：Gaussian他居然不支持谐振势二面角参数，只认周期势

没办法，只好硬着头皮把谐振势拟合成周期势。这个顶多拟合平衡位置附近的行为，不要指望能拟合的多贴切
(, 下载次数 Times of downloads: 33)
在转换RTP文件时，如果发现谐振势二面角，程序将使用最小二乘法拟合周期势，拟合范围为平衡位置±60°。会尝试一些常规的相位组合，取RMSE最低的一组：
[attach]113113[/attach]
(, 下载次数 Times of downloads: 31)
如果不满意，可以自定义相位重新拟合，但一般不怎么推荐，因为默认的相位拟合不动的一般要么有位阻要么有静电作用，强行拟合是不现实的，如果不太重要直接在平衡角度写个很大的数让它变成刚性的就得了。哎说起来这个功能应该很容易实现来着，下次更新加进去好了（）

再次回车，就如前述般直接生成gjf文件，同样可以直接进行计算！这里还会额外生成一个shell脚本，用于将拟合曲线的数据使用gnuplot绘图，检查拟合情况。

tips：对于需要旋转的二面角，去凑谐振势就完全是胡搞了。此时应当通过扫PES来进行拟合，笔者这里强烈安利ggdh老师的ztop，其torsionfit模块可以自动生成扫描所需的量子化学输入文件，还能以多种方式进行拟合，非常强大！

以上是主要功能，其他不太重要的玩法可以看readme或执行banetop -h打印帮助。
这个小玩具暂时写到这里了，后面有时间可能会再集成点oniom相关的功能- -。本软件是收费程序，售价5eV/人，熟人八折，自觉缴费。本程序的力场参数极大地依赖sobtop，因此如您使用本软件进行计算并用于发表，不要忘了按照sobtop主页的提示进行引用。有兴趣引本程序的话就引github吧，感觉这种玩具也没啥发表前途

作者
Author: 风起~ 时间: 2025-6-30 15:51
能不能支持一下ztop呀 ztop考虑的功能相对复杂一些
或者说可以把itp转为oniom格式的呀

作者
Author: wal 时间: 2025-7-1 13:12

风起~ 发表于 2025-6-30 15:51
能不能支持一下ztop呀 ztop考虑的功能相对复杂一些
或者说可以把itp转为oniom格式的呀

一开始是有考虑过支持ztop的，但是写着写着发现自己一个人精力不太够hh
ONIOM的话，如果我接下来有时间继续更新banetop，第一项任务肯定就是支持ONIOM。我个人目前还是用vmd与gview设置oniom输入文件的

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