|
本帖最后由 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` 配置文件即可。
- {
- "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;}"
- ],
- "energy": {
- "jobstep": "CCSD(T)-F12A",
- "property": "CCSD(T)-F12a"
- }
- }
- }
复制代码
上面是一个 `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 输出文件中有节点
- <jobstep command="CCSD(T)-F12A" commandset="CCSD">
复制代码
此处的 `command` 后的 `"CCSD(T)-F12A"` 即是需要填写到 `gmm.json` 的 `"jobstep"` 后的值。
再看该 `<jobstep>` 节点中,有节点
- <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. 单点
- %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 方法,似乎还可以,没有细测。
过渡态结构
能量变化
而 Molpro 直接算飞了……就放个最后的结构让大家笑笑吧。
五、总结
本文开发了一个联用 Gaussian、MOKIT 与 Molpro 的 Python 脚本,可以在 Molpro 的 Post-HF 计算级别下,结合 Gaussian 的 SCF 与优化算法,进行几何结构优化、过渡态搜索、频率计算与生成 IRC,并取得了比较好的效果。
六、附录
把代码的解析放在最后吧,良心级别的源码解读!
首先是解析 `gmm.json` 里写的内容,这里的变量已经在前面讲过了,直接看代码吧。
- with open("gmm.json") as fo:
- settings = json.load(fo)
- G_NPROC = settings["Gaussian"]["nproc"]
- G_MEMORY = settings["Gaussian"]["memory"]
- G_LEVEL = settings["Gaussian"]["level"]
- M_NPROC = settings["Molpro"]["nproc"]
- M_MEMORY = settings["Molpro"]["memory"]
- M_COMMANDS = settings["Molpro"]["commands"]
- M_ENERGY = settings["Molpro"]["energy"]
复制代码
这里为了让计算目录干净点,就在计算目录下新建一个 `tmp` 目录并进去,之后的计算都在这里完成。
- if not os.path.exists("tmp"):
- os.mkdir("tmp")
-
- os.chdir("tmp")
复制代码
接下来了构建 Gaussian 输入文件模板,以及单位转换,和一个元素周期表(太大了就不放了)。
- ROHF_INP = f"""%nproc={G_NPROC}
- %mem={G_MEMORY}
- %chk=mol_rohf.chk
- {G_LEVEL}
- ROHF TASK
- """
复制代码
接下来解析 Gaussian 喂进来的文件,读取原子数、电荷数、自旋多重度以及坐标,并填写到 Gaussian 输入文件模板中,写入到 `mol_rohf.gjf` 文件,开始 ROHF 计算。
- (layer, InputFile, OutputFile, MsgFile, FChkFile, MatElFile) = sys.argv[1:]
- atoms = []
- coords = []
- 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"
- with open("mol_rohf.gjf", "w") as fo:
- fo.write(ROHF_INP)
- print(">>> Starting ROHF Calculation...")
- os.system("g16 mol_rohf.gjf")
- os.system("formchk mol_rohf.chk")
- print(">>> ROHF Calculation Done!")
复制代码
ROHF 计算完成后,使用 MOKIT 中的 `fch2com` 生成 Molpro 的输入文件与轨道,并改个名字。再追加 `gmm.json` 里填写的 Molpro 计算命令。如果 Gaussian 要求能量的一阶梯度,脚本就会自动填写上 `force` 命令,就可以进行 Molpro 计算了。
- print(">>> Preparing Molpro input file...")
- os.system("fch2com mol_rohf.fchk")
- os.system("mv mol_rohf.com mol_post.com")
- # add addtional molpro calculation procedures
- for command in M_COMMANDS:
- os.system(f"echo '{command}' >> mol_post.com")
- if derivs == 1:
- os.system("echo '{force;}' >> mol_post.com")
- # call molpro to do post hartree fock calculation
- print(">>> Starting Molpro Calculation...")
- os.system(f"molpro -s -n {M_NPROC} -m {M_MEMORY} mol_post.com")
- print(">>> Molpro Calculation Done!")
复制代码
Molpro 计算完成后,就要解析 xml 文件了,解析过程其实就是按 `gmm.json` 中的填写的属性去定位,找出来还给 Gaussian 就行了。这块写得其实不是很满意,写 xpath 匹配会更优雅,但既然 work 了,就不改了。
- print(">>> Extracting Molpro output file ...")
- DomTree = parse("mol_post.xml")
- root = DomTree.documentElement
- jobstep_list = root.getElementsByTagName("jobstep")
- for jobstep in jobstep_list:
- if jobstep.getAttribute("command") == M_ENERGY["jobstep"]:
- for property in jobstep.getElementsByTagName("property"):
- if property.getAttribute("name") == "total energy" and property.getAttribute("method") == M_ENERGY["property"]:
- energy = float(property.getAttribute("value"))
- break
- break
- if derivs == 1:
- for jobstep in jobstep_list:
- if jobstep.getAttribute("command") == "FORCE":
- grad_text = jobstep.getElementsByTagName("gradient")[0].childNodes[0].data
- grad_lines = grad_text.split("\n")
- gradients = []
- for line in grad_lines:
- atom_grad = [float(x) for x in line.split()]
- if len(atom_grad) == 3:
- gradients.append(atom_grad)
- with open(OutputFile, "w") as fo:
- fo.writelines(f"{energy:20.12E}{0:20.12E}{0:20.12E}{0:20.12E}\n")
- if derivs == 1:
- for atom_grad in gradients:
- fo.writelines("".join([f"{g:20.12E}" for g in atom_grad]) + "\n")
复制代码
|
-
-
gmm.zip
218.57 KB, 下载次数 Times of downloads: 23
评分 Rate
-
查看全部评分 View all ratings
|