计算化学公社
标题: Gaussian-MOKIT-Molpro 联用方案 gmm.py [打印本页]
作者Author: mizu-bai 时间: 2022-11-22 22:06
标题: Gaussian-MOKIT-Molpro 联用方案 gmm.py
本帖最后由 mizu-bai 于 2024-1-23 14:26 编辑
一、前言
自从看了社长写的 Gaussian 与 xTB/ORCA 联用方案后，一直想把组里常用的 Molpro 与 Gaussian 接起来使用。苦于 Molpro 输入文件的模板过于难写，一直没有动手。后来了解了邹神的 MOKIT，里面的 fch2com 小工具可以方便地将 Gaussian 做 ROHF 计算的 fchk 文件转成 Molpro 的输入文件，同时也将基组和轨道传入，让 SCF 收敛简单了不少。笔者之前用 Molpro 做优化时，就常常遇到这个问题，尝试了挺多方法，但有时候实在还是 RHF 无法收敛，而几何结构优化遇到问题也已经习惯了，不收敛，多几个虚频都是常事。
由于我不怎么会写 Shell 脚本，最熟悉的是 Python 和 Perl，又考虑到了 MOKIT 推荐用 Intel oneAPI 套件编译，而里面又自带了个 Python，最终花了两三天时间，写了个 Python 脚本把 Gaussian，MOKIT 和 Molpro 这三个软件捏起来了。已在 GitHub 上以 BSD-2 协议开源 https://github.com/mizu-bai/Gaussian-MOKIT-Molpro，README.md 即英文版手册。
至于 Gaussian、MOKIT 和 Molpro 的安装就不多赘述了，如果像本文中只使用 MOKIT 中 `fch2inp` 的功能的话，就只编译 MOKIT 本体，不用去装 PySCF 那些了。
二、计算流程
简单介绍一下整体的计算流程。
1. 首先有一个通过 external 关键词调用外部脚本的 Gaussian 输入文件，这个文件控制执行的任务和分子结构等，注意这个脚本内要指定 `%nproc=1`；
2. 接下来 Gaussian 会调用笔者写的这个 Python 脚本，脚本会读取传入的文件里的分子的电荷、自旋多重度和坐标，并生成一个 ROHF 的计算文件，再调用另一个 Gaussian 去做 ROHF 计算；
3. ROHF 计算完成后，调用 MOKIT 中的 `fch2com` 转成 Molpro 输入文件和波函数文件，即用 Gaussian 先帮 Molpro 做了 ROHF 计算，再让 Molpro 去做 Post-HF 计算；
4. 脚本将用户自定义的计算命令追加到 fch2com 生成的 Molpro 计算文件中，并调用 Molpro 计算；
5. Molpro 计算完成后，会将能量、偶极矩与梯度保存在 csv 文件中，解析后返回给 Gaussian。
由于 Molpro 只支持闭壳层的 HF 和 MCSCF 的解析频率，再加上太懒了，不想解析 Molpro 的 `frequencies` 输出结果，干脆直接在 Gaussian 里 `freq=(num)` 好了。以及事实上 Molpro 也没多少方法支持解析梯度，笔者虽然解析了 Molpro 的 `force` 命令输出结果，Gaussian 以为这方法能算解析梯度，实际上梯度可能是 Molpro 数值算出来的。
三、使用说明
由于之前使用社长的 Gaussian 与 xTB/ORCA 联用的脚本时，计算级别都写在对应的 `*.sh` 文件中，每次使用都复制一份到当前目录下，还要按需修改。这次写脚本的时候就考虑把指定计算级别的部分拆分开来，于是就有了 `gmm.py` 和 `gmm.json` 两个文件，`gmm.py` 的内容不需要改动，所以可以放在一个公用可以访问到的地方，而计算目录下只要放一个 `gmm.json` 配置文件即可。
- {
- "Gaussian": {
- "nproc": 4,
- "memory": "4GB",
- "level": "#p rohf aug-cc-pvdz nosymm int=(nobasistransform)"
- },
- "Molpro": {
- "nproc": 8,
- "memory": "200M",
- "commands": [
- "basis=avdz",
- "{rhf;}",
- "{ccsd(t)-f12a;}"
- ]
- }
- }
上面是一个 `gmm.json` 配置文件的例子，是一个 json 文件，哪怕不清楚 json 文件的格式也应该能知道哪一块是干什么的。`Gaussian` 下面的 `nproc` 和 `memory` 指定了计算资源，而 `level` 指定了计算级别，注意一定要有 `nosymm` 和 `int=(nobasistransform)`，为什么要加这两个可以去 MOKIT 主页了解。`Molpro` 下面同样指定了计算资源，`commands` 里是需要追加到 `fch2com` 生成的输入文件里的计算任务，由于自定义基组不能匹配到辅助基组，所以在进行某些任务的时候还必须重新指定基组，做一次 HF 计算，再做 Post-HF 计算。当然这些 HF 计算基本都一圈收敛，耗费的时候可以忽略。如果调用 gmm.py 的 Gaussian 输入文件要求计算梯度，如 opt 或 force 等任务，则 gmm.py 会自动追加梯度计算命令于 Molpro 输入文件中。之后便可以在 Gaussian 输入文件中以 `external="python3 -u /path/to/gmm.py"` 运行任务，`gmm.py` 会自动读取计算目录下的 `gmm.json` 文件。
下面给几种计算任务输入文件的例子，另外附件里也会有一些开发过程中测试文件，有兴趣的读者可以自己算一下。
1. 单点
- %nproc=1
- %chk=H2.chk
- # external="python3 -u /path/to/gmm.py"
- H2
- 0 1
- H 0.00000000 0.00000000 -0.30000000
- H 0.00000000 0.00000000 0.30000000
2. 优化 + 频率
- %nproc=1
- %chk=H2.chk
- # opt=(nomicro) external="python3 -u /path/to/gmm.py"
- H2
- 0 1
- H 0.00000000 0.00000000 -0.30000000
- H 0.00000000 0.00000000 0.30000000
- --link1--
- %nproc=1
- %chk=H2.chk
- # freq=(num) geom=(allcheck) external="python3 -u /path/to/gmm.py"
3. 过渡态 + 频率，需要先做一次频率计算，再读取力常数优化过渡态，最后计算频率，相当于 rcfc
- %nproc=1
- %chk=H3.chk
- # freq=(num) external="python3 -u /path/to/gmm.py"
- H3
- 0 2
- H 0.00000000 0.00000000 1.33333333
- H 0.00000000 0.00000000 0.33333333
- H 0.00000000 0.00000000 -1.66666667
- --link1--
- %nproc=1
- %chk=H3.chk
- # opt=(nomicro,rcfc,ts,noeigen) geom=(allcheck) external="python3 -u /path/to/gmm.py"
- --link1--
- %nproc=1
- %chk=H3.chk
- # freq=(num) geom=(allcheck) external="python3 -u /path/to/gmm.py"
-
4. IRC，基于过渡态搜索完成后的频率计算结果
- %nproc=1
- %chk=H3.chk
- # irc=(rcfc,gradientonly) geom=(allcheck) external="python3 -u /path/to/gmm.py"
四、测试效果
笔者使用了社长的 Gaussian 与 ORCA 联用中 HCN 与 CNH 异构化反应的输入文件中的结构，分别使用脚本联用和单纯使用 Molpro 两种方式，尝试计算过渡态和 IRC，对比两者的优化算法差异。计算级别统一为 CCSD(T)-F12a/AVTZ，测试文件可以在附件里获取。
脚本联用时只在最初计算了力常数，而 Molpro 的默认设置则自动选择了 QSD 方法，每步计算一次 Hessian，但最终的结果令人大跌眼镜。脚本联用时，3 步能量基本稳定，7 步时最终收敛，虚频 1195.74 cm^-1，IRC 也很完美，社长在过渡态和 IRC 那篇文章中提到过只用能量梯度的 IRC 算法效果并不怎么样，但在 Gaussian 16 中默认选项换为了 EluerPC 方法，似乎还可以，没有细测。
(, 下载次数 Times of downloads: 19)
(, 下载次数 Times of downloads: 14)
(, 下载次数 Times of downloads: 15)
(, 下载次数 Times of downloads: 12)
而 Molpro 直接算飞了……就放个最后的结构让大家笑笑吧。
(, 下载次数 Times of downloads: 13)
五、总结
本文开发了一个联用 Gaussian、MOKIT 与 Molpro 的 Python 脚本，可以在 Molpro 的 Post-HF 计算级别下，结合 Gaussian 的 SCF 与优化算法，进行几何结构优化、过渡态搜索、频率计算与生成 IRC，并取得了比较好的效果。
作者Author: 鬼隐 时间: 2022-11-22 22:09
mizu好厉害！ 顶！这个package好实用！
贴贴
作者Author: mizu-bai 时间: 2022-11-26 14:32
Update
修改了 gmm.py 中追加 Molpro Post-HF 计算命令的代码，将每次计算的能量和梯度信息保存到 energy.csv 和 grad.csv 中，不用再在 gmm.json 指定 xml 需要解析的节点的属性值了，帖子中的 gmm.json 应修改成下面样子：
- {
- "Gaussian": {
- "nproc": 4,
- "memory": "4GB",
- "level": "#p rohf aug-cc-pvdz nosymm int=(nobasistransform)"
- },
- "Molpro": {
- "nproc": 8,
- "memory": "200M",
- "commands": [
- "basis=avdz",
- "{rhf;}",
- "{ccsd(t)-f12a;}"
- ]
- }
- }
最新版 gmm.py 代码在此
(, 下载次数 Times of downloads: 15)
作者Author: shenzp 时间: 2023-2-11 23:39
现在这个脚本好像不能把gaussian里的自定义基组放进生成的gjf文件中，请问怎么改一下能用自定义基组？
作者Author: mizu-bai 时间: 2023-2-12 14:08
本帖最后由 mizu-bai 于 2023-2-12 14:10 编辑
注意看源码中生成 Gaussian ROHF 输入文件那段
- with open(InputFile, "r") as fi:
- (atoms, derivs, charge, spin) = [int(x) for x in fi.readline().split()]
- ROHF_INP += f"{charge} {spin}\n"
- for i in range(0, atoms):
- arr = fi.readline().split()
- atom = ELEMENTS[int(arr[0])]
- coord = [f"{(float(x) / ANG2BOHR)}" for x in arr[1:4]]
- ROHF_INP += f" {atom}"
- ROHF_INP += " " * 4
- ROHF_INP += (" " * 4).join(coord)
- ROHF_INP += "\n"
- ROHF_INP += "\n"
在最后一个 `ROHF_INP += "\n"` 后添加一行
- ROHF_INP += "@basis_file.gbs"
注意这一行的缩进一定要和最后一个 `ROHF_INP += "\n"` 对齐，均为四个空格，并把自己的基组文件放到当前目录下。以及记得把 `gmm.json` 配置文件中的 Gaussian 计算级别中基组的部分换成 `gen`。
作者Author: shenzp 时间: 2023-2-12 19:12
好的，非常感谢！
作者Author: mizu-bai 时间: 2023-5-4 10:50
补一下 Benchmark 的数据
(, 下载次数 Times of downloads: 14)
(, 下载次数 Times of downloads: 14)
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |