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 这三个软件捏起来了。

至于 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 计算完成后，脚本解析 Molpro 输出的 xml 文件，返回能量和梯度给 Gaussian。

由于 Molpro 除了闭壳层的 HF 和 MCSCF 支持解析频率，再加上太懒了，不想解析 Molpro 的 `frequencies` 输出结果，干脆直接在 Gaussian 里 `freq=(num)` 好了。比较奇怪的一点是，通过本脚本让 Gaussian 基于能量梯度算频率居然比 Molpro 直接用算频率（也是数值的）要更快。以及事实上 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.         "energy": {
    
16.             "jobstep": "CCSD(T)-F12A",
    
17.             "property": "CCSD(T)-F12a"
    
18.         }
    
19.     }
    
20. }

复制代码

上面是一个 `gmm.json` 配置文件的例子，是一个 json 文件，哪怕不清楚 json 文件的格式也应该能知道哪一块是干什么的。`Gaussian` 下面的 `nproc` 和 `memory` 指定了计算资源，而 `level` 指定了计算级别，注意一定要有 `nosymm` 和 `int=(nobasistransform)`，为什么要加这两个可以去 MOKIT 主页了解。`Molpro` 下面同样指定了计算资源，`commands` 里是需要追加到 `fch2com` 生成的输入文件里的计算任务，由于自定义基组不能匹配到辅助基组，所以在进行某些任务的时候还必须重新指定基组，做一次 HF 计算，再做 Post-HF 计算。当然这些 HF 计算基本都一圈收敛，耗费的时候可以忽略。而 `energy` 段是最需要注意的，因为本脚本是通过匹配 Molpro 输出的 xml 文件获得需要的能量，所以需要确切的 xml 节点属性，笔者在下面列了一些常见的 Post-HF 计算需要匹配的内容，按需填写就好。如果没有的话，就先找个很简单的体系做对应级别的计算，先去找一下 xml 输出文件里的 `command` 属性为你的计算级别的 `<jobstep>` 节点，把 `command` 属性对应的字符串填到 `gmm.json` 配置文件中，再在该 `<jobstep>` 节点下找 `name` 属性为 `"total energy"` 的 `<property>` 节点，把 `method` 属性对应的字符串填写到 `gmm.json` 中即可。

下面以 CCSD(T)-F12a 级别为例，在 xml 输出文件中有节点

1. <jobstep command="CCSD(T)-F12A" commandset="CCSD">

复制代码

此处的 `command` 后的 `"CCSD(T)-F12A"` 即是需要填写到 `gmm.json` 的 `"jobstep"` 后的值。

再看该 `<jobstep>` 节点中，有节点

1. <property name="total energy" method="CCSD(T)-F12a" principal="true" stateSymmetry="1" stateNumber="1" value="-139.601945640719"/>

复制代码

此处的 `name` 属性为 `"total energy"`，便是我们需要的能量值，将后面的 `"CCSD(T)-F12a"` 填入 `gmm.json` 的 `"property"` 后即可。

之后便可以在 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，并取得了比较好的效果。

六、附录

把代码的解析放在最后吧，良心级别的源码解读！

首先是解析 `gmm.json` 里写的内容，这里的变量已经在前面讲过了，直接看代码吧。

1. with open("gmm.json") as fo:
   
2.     settings = json.load(fo)
   
3. 
   
4. G_NPROC    = settings["Gaussian"]["nproc"]
   
5. G_MEMORY   = settings["Gaussian"]["memory"]
   
6. G_LEVEL    = settings["Gaussian"]["level"]
   
7. 
   
8. M_NPROC    = settings["Molpro"]["nproc"]
   
9. M_MEMORY   = settings["Molpro"]["memory"]
   
10. M_COMMANDS = settings["Molpro"]["commands"]
    
11. M_ENERGY   = settings["Molpro"]["energy"]

复制代码

这里为了让计算目录干净点，就在计算目录下新建一个 `tmp` 目录并进去，之后的计算都在这里完成。

1. if not os.path.exists("tmp"):
   
2.     os.mkdir("tmp")
   
3.         
   
4. os.chdir("tmp")

复制代码

接下来了构建 Gaussian 输入文件模板，以及单位转换，和一个元素周期表（太大了就不放了）。

1. ROHF_INP = f"""%nproc={G_NPROC}
   
2. %mem={G_MEMORY}
   
3. %chk=mol_rohf.chk
   
4. {G_LEVEL}
   
5. 
   
6. ROHF TASK
   
7. 
   
8. """

复制代码

接下来解析 Gaussian 喂进来的文件，读取原子数、电荷数、自旋多重度以及坐标，并填写到 Gaussian 输入文件模板中，写入到 `mol_rohf.gjf` 文件，开始 ROHF 计算。

1. (layer, InputFile, OutputFile, MsgFile, FChkFile, MatElFile) = sys.argv[1:]
   
2. 
   
3. atoms = []
   
4. coords = []
   
5. with open(InputFile, "r") as fi:
   
6.     (atoms, derivs, charge, spin) = [int(x) for x in fi.readline().split()]
   
7.     ROHF_INP += f"{charge} {spin}\n"
   
8.     for i in range(0, atoms):
   
9.         arr = fi.readline().split()
   
10.         atom = ELEMENTS[int(arr[0])]
    
11.         coord = [f"{(float(x) / ANG2BOHR)}" for x in arr[1:4]]
    
12.         ROHF_INP += f" {atom}"
    
13.         ROHF_INP += " " * 4
    
14.         ROHF_INP += (" " * 4).join(coord)
    
15.         ROHF_INP += "\n"
    
16. 
    
17.     ROHF_INP += "\n"
    
18. 
    
19. with open("mol_rohf.gjf", "w") as fo:
    
20.     fo.write(ROHF_INP)
    
21. 
    
22. print(">>> Starting ROHF Calculation...")
    
23. os.system("g16 mol_rohf.gjf")
    
24. os.system("formchk mol_rohf.chk")
    
25. print(">>> ROHF Calculation Done!")

复制代码

ROHF 计算完成后，使用 MOKIT 中的 `fch2com` 生成 Molpro 的输入文件与轨道，并改个名字。再追加 `gmm.json` 里填写的 Molpro 计算命令。如果 Gaussian 要求能量的一阶梯度，脚本就会自动填写上 `force` 命令，就可以进行 Molpro 计算了。

1. print(">>> Preparing Molpro input file...")
   
2. os.system("fch2com mol_rohf.fchk")
   
3. os.system("mv mol_rohf.com mol_post.com")
   
4. 
   
5. # add addtional molpro calculation procedures
   
6. for command in M_COMMANDS:
   
7.     os.system(f"echo '{command}' >> mol_post.com")
   
8. 
   
9. if derivs == 1:
   
10.     os.system("echo '{force;}' >> mol_post.com")
    
11. 
    
12. # call molpro to do post hartree fock calculation
    
13. print(">>> Starting Molpro Calculation...")
    
14. os.system(f"molpro -s -n {M_NPROC} -m {M_MEMORY} mol_post.com")
    
15. print(">>> Molpro Calculation Done!")

复制代码

Molpro 计算完成后，就要解析 xml 文件了，解析过程其实就是按 `gmm.json` 中的填写的属性去定位，找出来还给 Gaussian 就行了。这块写得其实不是很满意，写 xpath 匹配会更优雅，但既然 work 了，就不改了。

1. print(">>> Extracting Molpro output file ...")
   
2. DomTree = parse("mol_post.xml")
   
3. root = DomTree.documentElement
   
4. jobstep_list = root.getElementsByTagName("jobstep")
   
5. 
   
6. for jobstep in jobstep_list:
   
7.     if jobstep.getAttribute("command") == M_ENERGY["jobstep"]:
   
8.         for property in jobstep.getElementsByTagName("property"):
   
9.             if property.getAttribute("name") == "total energy" and property.getAttribute("method") == M_ENERGY["property"]:
   
10.                 energy = float(property.getAttribute("value"))
    
11.                 break
    
12. 
    
13.         break
    
14. 
    
15. if derivs == 1:
    
16.     for jobstep in jobstep_list:
    
17.         if jobstep.getAttribute("command") == "FORCE":
    
18.             grad_text = jobstep.getElementsByTagName("gradient")[0].childNodes[0].data
    
19.             grad_lines = grad_text.split("\n")
    
20.             gradients = []
    
21.             for line in grad_lines:
    
22.                 atom_grad = [float(x) for x in line.split()]
    
23.                 if len(atom_grad) == 3:
    
24.                     gradients.append(atom_grad)
    
25. 
    
26. with open(OutputFile, "w") as fo:
    
27.     fo.writelines(f"{energy:20.12E}{0:20.12E}{0:20.12E}{0:20.12E}\n")
    
28.     if derivs == 1:
    
29.     for atom_grad in gradients:
    
30.         fo.writelines("".join([f"{g:20.12E}" for g in atom_grad]) + "\n")
    

复制代码

鬼隐

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: 7)

2022-11-26 14:31 上传 Uploaded

点击下载
Click to download