计算化学公社

 找回密码 Forget password
 注册 Register
Views: 790|回复 Reply: 2

[辅助/分析程序] Gaussian-MOKIT-Molpro 联用方案 gmm.py

[复制链接 Copy URL]

18

帖子

2

威望

578

eV
积分
636

Level 4 (黑子)

发表于 Post on 2022-11-22 22:06:22 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 mizu-bai 于 2022-11-22 22:05 编辑

一、前言

自从看了社长写的 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. H2

  5. 0 1
  6. H        0.00000000    0.00000000   -0.30000000
  7. 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. H2

  5. 0 1
  6. H        0.00000000    0.00000000   -0.30000000
  7. H        0.00000000    0.00000000    0.30000000


  8. --link1--
  9. %nproc=1
  10. %chk=H2.chk
  11. # freq=(num) geom=(allcheck) external="python3 -u /path/to/gmm.py"

复制代码

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

  1. %nproc=1
  2. %chk=H3.chk
  3. # freq=(num) external="python3 -u /path/to/gmm.py"

  4. H3

  5. 0 2
  6. H    0.00000000    0.00000000    1.33333333
  7. H    0.00000000    0.00000000    0.33333333
  8. H    0.00000000    0.00000000   -1.66666667


  9. --link1--
  10. %nproc=1
  11. %chk=H3.chk
  12. # opt=(nomicro,rcfc,ts,noeigen) geom=(allcheck) external="python3 -u /path/to/gmm.py"


  13. --link1--
  14. %nproc=1
  15. %chk=H3.chk
  16. # freq=(num) geom=(allcheck) external="python3 -u /path/to/gmm.py"

复制代码

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 方法,似乎还可以,没有细测。

过渡态结构

过渡态结构


能量变化

能量变化


HCN-Vib.png


HCN-IRC.png


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

Molpro-Fail.png


五、总结

本文开发了一个联用 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. G_NPROC    = settings["Gaussian"]["nproc"]
  4. G_MEMORY   = settings["Gaussian"]["memory"]
  5. G_LEVEL    = settings["Gaussian"]["level"]

  6. M_NPROC    = settings["Molpro"]["nproc"]
  7. M_MEMORY   = settings["Molpro"]["memory"]
  8. M_COMMANDS = settings["Molpro"]["commands"]
  9. 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. ROHF TASK

  6. """
复制代码

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

  1. (layer, InputFile, OutputFile, MsgFile, FChkFile, MatElFile) = sys.argv[1:]

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

  15.     ROHF_INP += "\n"

  16. with open("mol_rohf.gjf", "w") as fo:
  17.     fo.write(ROHF_INP)

  18. print(">>> Starting ROHF Calculation...")
  19. os.system("g16 mol_rohf.gjf")
  20. os.system("formchk mol_rohf.chk")
  21. 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. # add addtional molpro calculation procedures
  5. for command in M_COMMANDS:
  6.     os.system(f"echo '{command}' >> mol_post.com")

  7. if derivs == 1:
  8.     os.system("echo '{force;}' >> mol_post.com")

  9. # call molpro to do post hartree fock calculation
  10. print(">>> Starting Molpro Calculation...")
  11. os.system(f"molpro -s -n {M_NPROC} -m {M_MEMORY} mol_post.com")
  12. 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. for jobstep in jobstep_list:
  6.     if jobstep.getAttribute("command") == M_ENERGY["jobstep"]:
  7.         for property in jobstep.getElementsByTagName("property"):
  8.             if property.getAttribute("name") == "total energy" and property.getAttribute("method") == M_ENERGY["property"]:
  9.                 energy = float(property.getAttribute("value"))
  10.                 break

  11.         break

  12. if derivs == 1:
  13.     for jobstep in jobstep_list:
  14.         if jobstep.getAttribute("command") == "FORCE":
  15.             grad_text = jobstep.getElementsByTagName("gradient")[0].childNodes[0].data
  16.             grad_lines = grad_text.split("\n")
  17.             gradients = []
  18.             for line in grad_lines:
  19.                 atom_grad = [float(x) for x in line.split()]
  20.                 if len(atom_grad) == 3:
  21.                     gradients.append(atom_grad)

  22. with open(OutputFile, "w") as fo:
  23.     fo.writelines(f"{energy:20.12E}{0:20.12E}{0:20.12E}{0:20.12E}\n")
  24.     if derivs == 1:
  25.     for atom_grad in gradients:
  26.         fo.writelines("".join([f"{g:20.12E}" for g in atom_grad]) + "\n")
复制代码

gmm.zip

218.57 KB, 下载次数 Times of downloads: 23

评分 Rate

参与人数
Participants 12
威望 +1 eV +54 收起 理由
Reason
ggsddu + 4 好物!
ChemG + 5 GJ!
zjxitcc + 5 精品内容
dantevinsky + 5 GJ!
丁越 + 5 赞!
冰释之川 + 5 Gjj的男朋友好腻害
sobereva + 1
元气蛋 + 5 好物!
卡开发发 + 5 хорошо!
hebrewsnabla + 5 精品内容
鬼隐 + 5 好物!
biogon + 5 GJ!

查看全部评分 View all ratings

20

帖子

0

威望

1046

eV
积分
1066

Level 4 (黑子)

发表于 Post on 2022-11-22 22:09:02 | 显示全部楼层 Show all
mizu好厉害! 顶!这个package好实用!
贴贴

18

帖子

2

威望

578

eV
积分
636

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2022-11-26 14:32:04 | 显示全部楼层 Show all
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)

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 03:22 , Processed in 0.334612 second(s), 26 queries .

快速回复 返回顶部 返回列表 Return to list