Gaussian-MOKIT-Molpro 联用方案 gmm.py

mizu-bai

一、前言

自从看了社长写的 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` 配置文件即可。

1. {
   
2.     "Gaussian": {
   
3.         "nproc": 4,
   
4.         "memory": "4GB",
   
5.         "level": "#p rohf aug-cc-pvdz nosymm int=(nobasistransform)"
   
6.     },
   
7.     "Molpro": {
   
8.         "nproc": 8,
   
9.         "memory": "200M",
   
10.         "commands": [
    
11.             "basis=avdz",
    
12.             "{rhf;}",
    
13.             "{ccsd(t)-f12a;}"
    
14.         ]
    
15.     }
    
16. }

复制代码

上面是一个 `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. 单点

1. %nproc=1
   
2. %chk=H2.chk
   
3. # external="python3 -u /path/to/gmm.py"
   
4. 
   
5. H2
   
6. 
   
7. 0 1
   
8. H        0.00000000    0.00000000   -0.30000000
   
9. H        0.00000000    0.00000000    0.30000000

复制代码

2. 优化 + 频率

1. %nproc=1
   
2. %chk=H2.chk
   
3. # opt=(nomicro) external="python3 -u /path/to/gmm.py"
   
4. 
   
5. H2
   
6. 
   
7. 0 1
   
8. H        0.00000000    0.00000000   -0.30000000
   
9. H        0.00000000    0.00000000    0.30000000
   
10. 
    
11. 
    
12. --link1--
    
13. %nproc=1
    
14. %chk=H2.chk
    
15. # freq=(num) geom=(allcheck) external="python3 -u /path/to/gmm.py"
    
16. 
    

复制代码

3. 过渡态 + 频率，需要先做一次频率计算，再读取力常数优化过渡态，最后计算频率，相当于 rcfc

1. %nproc=1
   
2. %chk=H3.chk
   
3. # freq=(num) external="python3 -u /path/to/gmm.py"
   
4. 
   
5. H3
   
6. 
   
7. 0 2
   
8. H    0.00000000    0.00000000    1.33333333
   
9. H    0.00000000    0.00000000    0.33333333
   
10. H    0.00000000    0.00000000   -1.66666667
    
11. 
    
12. 
    
13. --link1--
    
14. %nproc=1
    
15. %chk=H3.chk
    
16. # opt=(nomicro,rcfc,ts,noeigen) geom=(allcheck) external="python3 -u /path/to/gmm.py"
    
17. 
    
18. 
    
19. --link1--
    
20. %nproc=1
    
21. %chk=H3.chk
    
22. # freq=(num) geom=(allcheck) external="python3 -u /path/to/gmm.py"
    
23. 
    

复制代码

4. IRC，基于过渡态搜索完成后的频率计算结果

1. %nproc=1
   
2. %chk=H3.chk
   
3. # 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 方法，似乎还可以，没有细测。

而 Molpro 直接算飞了……就放个最后的结构让大家笑笑吧。

五、总结

本文开发了一个联用 Gaussian、MOKIT 与 Molpro 的 Python 脚本，可以在 Molpro 的 Post-HF 计算级别下，结合 Gaussian 的 SCF 与优化算法，进行几何结构优化、过渡态搜索、频率计算与生成 IRC，并取得了比较好的效果。

鬼隐

mizu好厉害！ 顶！这个package好实用！
贴贴

mizu-bai

Update

修改了 gmm.py 中追加 Molpro Post-HF 计算命令的代码，将每次计算的能量和梯度信息保存到 energy.csv 和 grad.csv 中，不用再在 gmm.json 指定 xml 需要解析的节点的属性值了，帖子中的 gmm.json 应修改成下面样子：

1. {
   
2.     "Gaussian": {
   
3.         "nproc": 4,
   
4.         "memory": "4GB",
   
5.         "level": "#p rohf aug-cc-pvdz nosymm int=(nobasistransform)"
   
6.     },
   
7.     "Molpro": {
   
8.         "nproc": 8,
   
9.         "memory": "200M",
   
10.         "commands": [
    
11.             "basis=avdz",
    
12.             "{rhf;}",
    
13.             "{ccsd(t)-f12a;}"
    
14.         ]
    
15.     }
    
16. }

复制代码

最新版 gmm.py 代码在此 gmm.py (3.98 KB, 下载次数 Times of downloads: 15)

2022-11-26 14:31 上传 Uploaded

点击下载
Click to download

shenzp

现在这个脚本好像不能把gaussian里的自定义基组放进生成的gjf文件中，请问怎么改一下能用自定义基组？

mizu-bai

shenzp 发表于 2023-2-11 23:39
现在这个脚本好像不能把gaussian里的自定义基组放进生成的gjf文件中，请问怎么改一下能用自定义基组？

注意看源码中生成 Gaussian ROHF 输入文件那段

1. 
   
2. with open(InputFile, "r") as fi:
   
3.     (atoms, derivs, charge, spin) = [int(x) for x in fi.readline().split()]
   
4.     ROHF_INP += f"{charge} {spin}\n"
   
5.     for i in range(0, atoms):
   
6.         arr = fi.readline().split()
   
7.         atom = ELEMENTS[int(arr[0])]
   
8.         coord = [f"{(float(x) / ANG2BOHR)}" for x in arr[1:4]]
   
9.         ROHF_INP += f" {atom}"
   
10.         ROHF_INP += " " * 4
    
11.         ROHF_INP += (" " * 4).join(coord)
    
12.         ROHF_INP += "\n"
    
13. 
    
14.     ROHF_INP += "\n"
    

复制代码

在最后一个 `ROHF_INP += "\n"` 后添加一行

1. 
   
2.     ROHF_INP += "@basis_file.gbs"
   

复制代码

注意这一行的缩进一定要和最后一个 `ROHF_INP += "\n"` 对齐，均为四个空格，并把自己的基组文件放到当前目录下。以及记得把 `gmm.json` 配置文件中的 Gaussian 计算级别中基组的部分换成 `gen`。

shenzp

mizu-bai 发表于 2023-2-12 14:08
注意看源码中生成 Gaussian ROHF 输入文件那段

好的，非常感谢！

mizu-bai

补一下 Benchmark 的数据