利用MOKIT从PySCF向其他量化程序传轨道

2022-12-07 15:22:13 浏览数 (1)

MOKIT是免费、开源的轨道转换和接口程序,提供各种小程序和模块,能够在常见量子化学软件间传递分子轨道。其中的automr程序可以进行多参考态方法的自动化、黑箱式计算,详细介绍见《自动做多参考态计算的程序MOKIT》。近期笔者和另一开发者wsr在MOKIT程序中加入了fchk(),py2molpro,py2molcas,py2qchem等模块,可用于从PySCF程序向其他量子化学程序传递分子轨道。尤其是通过fchk()产生.fch文件,可方便地用于轨道可视化、波函数分析。

PySCF作为一款免费、开源的量子化学程序,如今已有众多用户。本公众号上发表过许多相关的安装和使用教程,例如 《PySCF程序包平均场计算的一些收敛技巧》(由孙启明本人撰写) 《离线安装PySCF-1.7.6》 《离线安装PySCF-2.x》

若能较方便地将PySCF与其他量子化学程序联用,做复杂方法的计算就会更得心应手。当然,此处的“联用”不是指简单地复制坐标,还要实现基组数据和分子轨道系数的正确格式转换,直接生成目标程序的输入文件和轨道文件,让其在计算时可以自动读入轨道。举三个例子:

(1)在PySCF中先做完CASSCF计算,到下一个程序里进行CASPT2计算,我们希望在下一个程序中CASSCF能极速收敛(指1圈收敛,或几圈收敛但能量不再下降),迅速进入CASPT2环节,节约计算时间。甚至直接以CASCI替代CASSCF轨道优化。 (2)目标程序的HF/DFT难以收敛,或能收敛但不支持检验波函数稳定性,那么我们可以先用PySCF做完HF/DFT计算,再传轨道至目标程序。(为什么不用Gaussian算完了传轨道给其他程序:因为Gaussian是商业收费程序,有的课题组/机构没买) (3)自己基于PySCF开发新方法,无现有程序对应,但希望正确地传轨道至下一个量化程序进行后续计算。

话不多说,直接看一个水分子的简单示例就明白。

1. PySCF传轨道给BDF

代码语言:javascript复制
from pyscf import gto, scf
from py2bdf import py2bdf

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2bdf(mf, 'h2o.inp')

产生h2o.inp,h2o.scforb和H2O.BAS三个文件。

2. PySCF传轨道给Dalton

代码语言:javascript复制
from pyscf import gto, scf
from py2dalton import py2dalton

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2dalton(mf, 'h2o.dal')

产生h2o.dal和h2o.mol两个文件。

3. PySCF传轨道给GAMESS

代码语言:javascript复制
from pyscf import gto, scf
from py2gms import py2gms

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2gms(mf, 'h2o.inp')

产生GAMESS输入文件h2o.inp,包含坐标、基组和轨道信息。

4. PySCF传轨道给OpenMolcas

代码语言:javascript复制
from pyscf import gto, scf
from py2molcas import py2molcas

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2molcas(mf, 'h2o.input')

产生h2o.input和h2o.INPORB文件。

5. PySCF传轨道给Molpro

代码语言:javascript复制
from pyscf import gto, scf
from py2molpro import py2molpro

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2molpro(mf, 'h2o.com')

产生h2o.com和h2o.a文件(含Alpha轨道)。如果是UHF/UDFT类型的计算,还会产生h2o.b文件(含Beta轨道)。读取.a和.b的关键词已在h2o.com中写好。

6. PySCF传轨道给ORCA

代码语言:javascript复制
from pyscf import gto, scf
from py2orca import py2orca

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2orca(mf, 'h2o.inp')

产生h2o.inp和h2o.mkl文件。如果检测到当前系统上有orca_2mkl小程序,则会自动将h2o.mkl转化为h2o.gbw文件。

7. PySCF传轨道给PSI4

代码语言:javascript复制
from pyscf import gto, scf
from py2psi import py2psi

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2psi(mf, 'h2o.inp')

产生h2o.inp和h2o.A文件(含Alpha轨道)。如果是UHF/UDFT类型的计算,还会产生h2o.B文件(含Beta轨道)。读取.A和.B文件的关键词已在h2o.inp中写好。

8. PySCF传轨道给Q-Chem

代码语言:javascript复制
from pyscf import gto, scf
from py2qchem import py2qchem

mol = gto.M(atom='''
O  -0.49390246   0.93902438   0.0
H   0.46609754   0.93902438   0.0
H  -0.81435705   1.84396021   0.0
''',
basis='cc-pVDZ')

mf = scf.RHF(mol).run()
py2qchem(mf, 'h2o.in')

产生h2o.in文件和h2o文件夹(内含轨道文件)。若当前机器存在环境变量$QCSCRATCH(即Q-Chem约定的临时文件存放目录),则h2o文件夹会被自动移入$QCSCRATCH/目录下,Q-Chem做计算时会自动识别之。例如运行

代码语言:javascript复制
qchem h2o.in h2o.out h2o

即可发现Q-Chem的RHF计算2圈收敛,能量无法再降低。若不存在$QCSCRATCH,则h2o文件夹仍留在当前目录下。

9. PySCF传轨道给Gaussian

代码语言:javascript复制
from pyscf import gto
from py2fch_direct import fchk

mol = gto.M(atom='O 0.0 0.0 0.1; H 0.0 0.0 1.0',
        basis='cc-pvtz',
        charge=-1,
        cart=True
        ).build()
mf = mol.RHF().run()
fchk(mf, 'test.fch')

相当于直接产生高斯fch文件。注意到这个模块的名称与上述其他模块不同,这是因为PSI4程序里有个叫fchk()的模块能够在计算完后导出fch文件,因此我们沿用了这个模块名称,希望用过PSI4的人都能对这个名称感到熟悉。这里我们再举一个传递CASSCF轨道的简单示例,算一个三重态氧气的CASSCF(6e,8o)

代码语言:javascript复制
from pyscf import gto, mcscf
from py2fch_direct import fchk

mol = gto.M(
    atom = 'O 0 0 0; O 0 0 1.2',
    basis = 'ccpvdz',
    spin = 2)

mf = mol.RHF().run()
mc = mcscf.CASSCF(mf, 6, 8).run() #CAS(6o,8e)
fchk(mc, 'O2_cas6o8e.fch')

注意这里我们其实提前看过O2的ROHF轨道,或对O2的分子轨道十分熟悉才能直接写(6e,8o),实际上是经过了观看、挑选轨道的过程。上述操作传的是CASSCF轨道,而更常用的是CASSCF自然轨道。这需要把最后两行改成如下四行

代码语言:javascript复制
mc = mcscf.CASSCF(mf, 6, 8) #CAS(6o,8e)
mc.natorb = True
mc.kernel()
fchk(mc, 'O2_cas6o8e_NO.fch')

这样O2_cas6o8e.fch文件里就有了CASSCF自然轨道和轨道占据数。若还想将CASSCF密度写进fch文件,则需将上述最后一行改为

代码语言:javascript复制
fchk(mc, 'O2_cas6o8e_NO.fch', density=True)

当然,这两种轨道对应的CASSCF能量都是一样的,我们可以用Gaussian检验一下,gjf文件如下

代码语言:javascript复制
%chk=O2_cas6o8e.chk
%mem=2GB
%nprocshared=2
#p CASSCF(8,6) chkbasis nosymm int=nobasistransform guess=read geom=allcheck

将fch转化为chk文件

代码语言:javascript复制
unfchk O2_cas6o8e.fch O2_cas6o8e.chk

提交计算

代码语言:javascript复制
g16 O2_cas6o8e.gjf &

可以发现2圈收敛且能量没有变化。O2_cas6o8e_NO.fch对应的计算就不再展示了,操作类似。

其他事项

为阅读简便起见,上述示例几乎都以水分子的RHF/cc-pVDZ计算为例,实际使用时不局限于此。读者需在生成的输入文件末尾加上自己需要的关键词以进行后续计算。另外注意在计算中请不要开启对称性,以免传轨道时有数据不对应问题。

MOKIT程序可在

代码语言:javascript复制
https://gitlab.com/jxzou/mokit

下载。主页上可以点鼠标下载源代码压缩包,也可以在下载框的Previous Artifacts中选择相应的Linux预编译版,解压后写好环境变量即可使用。注意Windows预编译版不支持本文功能,内含的是Gaussian与其他量化程序传轨道的小程序。若网页上所列的Linux预编译版不满足您的需求,可下载源代码自行编译,亦较为简单,见网页上的README_zh.md说明。有使用问题或建议欢迎留言,或直接在程序issues区提个issue。

0 人点赞