本文介绍如何使用MOKIT从ORCA向其他量化程序传轨道,有以下可能的用途:
(1)在ORCA中进行了RIJK或RIJCOSX加速的大体系HF/DFT计算,想传轨道给其他程序进行后续计算,或想产生fch文件方便可视化。
(2)在ORCA中进行了CCSD等post-HF计算,并产生了(自旋)自然轨道,想产生fch文件方便可视化。
(3)在ORCA中进行了CASSCF计算,想传CASSCF轨道给其他程序进行后续计算(例如MC-PDFT)。
(4)有些复杂体系可能ORCA可以收敛出特殊的SCF解,而目标程序难以得到,可以传轨道给目标程序。
为简洁起见,本文以水分子的RHF/def2-TZVP计算为例,输入文件h2o.inp内容如下
代码语言:javascript复制%pal nprocs 2 end
%maxcore 2000
! RHF def2-TZVP VeryTightSCF RIJCOSX def2/J defgrid3
* xyz 0 1
O 0.74608908 -0.25872442 0.0
H 1.70608908 -0.25872442 0.0
H 0.42563449 0.64621141 0.0
*
VeryTightSCF
是笔者的个人喜好,但做实际计算至少应采用TightSCF
。defgrid3
表示使用尽可能高的COSX格点。该例使用了RIJCOSX加速,但这对轨道系数影响非常小,更何况用了比较精细的格点和严格的收敛限,预期得到的轨道与一个传统的RHF/def2-TZVP计算的轨道极其接近。用ORCA算完后运行orca_2mkl h2o -mkl
转化生成h2o.mkl
文件,下面我们以该文件为例展示如何传轨道,读者若使用本文功能,请务必顺带阅读文末注意事项。
1. ORCA传轨道给AMESP
代码语言:javascript复制mkl2amo h2o.mkl
会产生h2o.amo和h2o.aip两个文件,分别含基组和轨道数据,坐标及关键词。为了让AMESP读入轨道,需要运行
代码语言:javascript复制a2m h2o.amo
即将h2o.amo文本文件转化为二进制文件h2o.mo,其中a2m
是AMESP自带的小程序。
2. ORCA传轨道给BDF
代码语言:javascript复制mkl2bdf h2o.mkl
会产生H2O.BAS,h2o.inp和h2o.scforb三个文件,分别含基组数据,坐标及关键词,轨道数据。
3. ORCA传轨道给CFOUR
代码语言:javascript复制mkl2cfour h2o.mkl
会产生ZMAT, GENBAS和OLDMOS三个文件,注意CFOUR所用文件名是固定的。mkl2cfour目前有个缺点:在书写ORCA输入文件时需要使用CFOUR输出文件里体系的直角坐标,产生的mkl文件才能给mkl2cfour使用。这是因为CFOUR会旋转分子朝向,只有使用CFOUR输出文件里体系的直角坐标,其已经被旋转至程序标准朝向,才往往不会再被旋转了。如果遇到CFOUR的SCF不收敛,显然此时有输出文件,这点倒不成问题。
4. ORCA传轨道给Dalton
代码语言:javascript复制mkl2dal h2o.mkl
会产生h2o.dal和h2o.mol两个文件。
5. ORCA传轨道给GAMESS
代码语言:javascript复制mkl2inp h2o.mkl
会产生h2o.inp文件,内含坐标,基组数据和一些简单的关键词。
6. ORCA传轨道给Molpro
代码语言:javascript复制mkl2com h2o.mkl
会产生h2o.a和h2o.com文件,前者含Alpha轨道,后者含坐标,基组和关键词。如果是UHF计算,还会产生h2o.b文件(含Beta轨道)。
7. ORCA传轨道给OpenMolcas
代码语言:javascript复制mkl2inporb h2o.mkl
会产生h2o.INPORB和h2o.input文件。
8. ORCA传轨道给PSI4
代码语言:javascript复制mkl2psi h2o.mkl
会产生h2o.A和h2o.inp文件。前者含Alpha轨道,后者含坐标,基组和关键词。如果是UHF类型计算,还会有h2o.B文件(内含Beta轨道)。
9. ORCA传轨道给PySCF
代码语言:javascript复制mkl2py h2o.mkl
会产生h2o.fch和h2o.py文件。在运行.py文件时会从.fch文件里读取轨道。
10. ORCA传轨道给Q-Chem
代码语言:javascript复制mkl2qchem h2o.mkl
产生h2o.in
文件和一个h2o
文件夹。如果机器上有定义环境变量$QCSCRATCH
,则h2o文件夹会自动被移至Q-Chem临时文件目录下,屏幕上提示你可以运行
qchem h2o.in h2o.out h2o
提交Q-Chem任务,最后一个参数表示读取h2o/目录下的轨道文件。若未定义$QCSCRATCH
,h2o文件夹则放在当前目录下,读者需要时自行移动。
11. ORCA传轨道给Gaussian
该功能较重要,有几种不同使用方式,此处重点介绍。
11.1 使用mkl2gjf小程序
代码语言:javascript复制mkl2gjf h2o.mkl
会产生h2o.gjf文件,内含坐标、基组信息和轨道信息。注意Gaussian是支持从gjf文件中读取轨道的(关键词guess=cards
)。
11.2 使用mkl2fch小程序
11.2.1 适用于常见HF/DFT/CASSCF轨道
代码语言:javascript复制mkl2fch h2o.mkl
若检测到当前目录下无h2o.fch文件,则会从零创建h2o.fch;若已存在h2o.fch文件,则仅将轨道信息写入h2o.fch。
11.2.2 适用于UNO、MP2、CCSD、CASSCF等自然轨道
此处以RI-MP2为例,输入文件计算级别部分为
代码语言:javascript复制! RHF def2-TZVP VeryTightSCF RIJCOSX def2/J defgrid3 RI-MP2 def2-TZVP/C
%mp2
density unrelaxed
natorbs true
end
关键词unrelaxed
表示计算RI-MP2非弛豫密度对应的自然轨道,若改成relaxed
即计算弛豫密度对应的自然轨道,笔者更推荐使用前者进行波函数分析。正常计算完后有一个文件叫h2o.mp2nat,它其实是gbw格式,运行如下几步
mv h2o.mp2nat h2o_mp2no.gbw
orca_2mkl h2o_mp2no -mkl
mkl2fch h2o_mp2no.mkl h2o_mp2no.fch -no
即生成含MP2自然轨道和轨道占据数的h2o_mp2no.fch文件。参数-no
意为将自然轨道占据数写入生成的fch文件。用GaussView打开该文件可视化轨道时,箭头旁的数值不是轨道能量,而是轨道占据数;用Multiwfn打开该文件可视化轨道时,Ene处数值为轨道占据数,Occ处数值无意义。
11.2.3 适用于UMP2、UCCSD(自旋)自然轨道
这类计算输入文件较复杂,往往涉及到对称破缺单重态和检验波函数稳定性等问题,计算实例见笔者在计算化学公社论坛上的回帖
代码语言:javascript复制http://bbs.keinsci.com/thread-29762-1-1.html
这里我们假设已经获得了.mdci.nat文件,接着执行
代码语言:javascript复制mv O_singlet_o.mdci.nat O_singlet_CCSD.gbw # 实际上是个gbw文件
orca_2mkl O_singlet_CCSD -mkl # 产生mkl
mkl2fch O_singlet_CCSD.mkl O_singlet_CCSD.fch -nso # 产生fch文件
参数-nso
意为将alpha、beta自旋自然轨道占据数写入生成的fch文件。打开该文件可视化轨道时轨道占据数含义同11.2.2。自旋自然轨道有alpha、beta两套,有的读者可能认为不便分析,但自然轨道只有1套。为获得UCCSD自然轨道,可启动Python,运行
from mokit.lib.lo import gen_no_from_nso
gen_no_from_nso(fchname='O_singlet_CCSD.fch')
产生UCCSD自然轨道文件O_singlet_CCSD_NO.fch,只含一套轨道。 《在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例》http://sobereva.com/403
注意事项
- 目前ORCA中开启点群对称性做计算只能旋转分子朝向、识别点群和不可约表示,并不能利用对称性加速计算,所以很少使用。是否开启点群对称性做计算不影响本文功能。但传轨道后在目标程序中应关闭对称性(生成的输入文件里写了这些关键词),笔者没有让各个小程序支持点群对称性的意向。
- 以上各个小程序都会产生目标程序输入文件,内含坐标和基组数据,不需要用户再去手动书写基组名称。强烈推荐用户使用该文件进行计算,既能免去手写基组的麻烦,也能保证传轨道时对应性更好。
- ORCA仅支持球谐型基函数,不支持笛卡尔型基函数,因此传轨道后产生的文件采用的也是球谐型基函数进行计算。
- 若ORCA的SCF有上千a.u.的剧烈振荡,很可能是出现了基函数线性依赖导致的,此时即使侥幸收敛了能量也未必靠谱,需要在输入文件里加上
%scf
sthresh 1e-6
end
此阈值是Gaussian和GAMESS处理基函数线性依赖的默认值,比较靠谱。
- 对于ORCA产生的.mp2nat,.mdci.nat和.nbo这些实际上是gbw格式的文件,转化成mkl文件前除了修改文件名外,还有一种办法
orca_2mkl h2o.mp2nat h2o_mp2no.mkl -mkl -anyorbs
就是两个文件名参数都带上扩展名,然后在最后加-anyorbs
参数,可省去改名麻烦。
- mkl文件有个缺点是不含赝势信息,笔者曾在ORCA官方论坛上发帖建议未来ORCA在mkl文件中增加赝势信息,但ORCA开发者表示目前没兴趣。若读者在计算中使用全电子基组,自然无此问题;若用了赝势,按上文操作产生其他量化程序的文件不会含赝势信息,即使轨道系数正确,SCF也会剧烈振荡。这里笔者推荐一种解决办法:
Step 1. 先用Gaussian算个几秒的任务 这里以二甲基锌(ZnMe2)为例,假设我们对Zn用SDD,gjf内容如下
代码语言:javascript复制%chk=ZnMe2.chk
#p RHF genecp nosymm int=nobasistransform guess(only,save)
title
0 1
[ZnMe2坐标,略]
H C 0
cc-pVDZ
****
Zn 0
SDD
****
Zn 0
SDD
关键词only
表示构建SCF初猜后即退出程序,即使大体系耗时也只有秒级别。若读者想让Gaussian正常算完,获得收敛的RHF轨道,可以删去guess(only,save)
。正常结束后获得ZnMe2.chk
文件。
Step 2. 获得含赝势的fch文件和ORCA相关文件 执行
代码语言:javascript复制fch2mkl ZnMe2.chk
随即产生ZnMe2.fch,ZnMe2_o.inp和ZnMe2_o.mkl文件,现在我们就有了一个含赝势信息的fch文件。ZnMe2_o.mkl含有来自Gaussian的轨道,除非有特殊计算目的,否则一般推荐使用该mkl文件,即运行
代码语言:javascript复制orca_2mkl ZnMe2_o -gbw
将mkl转化为gbw文件。我们用ZnMe2_o.inp文件在ORCA中做目标计算,最后把轨道传回已存在的fch文件
代码语言:javascript复制orca_2mkl ZnMe2_o -mkl # 假设ZnMe2_o.gbw已被更新,重新产生mkl
cp ZnMe2.fch ZnMe2_o.fch # 备份一下,以防覆盖
mkl2fch ZnMe2_o.mkl # 将轨道传回ZnMe2_o.fch
后续可以使用fch2inp
,fch2inporb
,fch2com
,bas_fch2py
等小程序传给其他量化程序做计算,且均含赝势信息。mkl2fch
检测到有fch文件存在时会实直接使用该文件,而非从零生成,因此避免了赝势缺失问题。注意不能跳过Step 1随便塞一个fch文件给mkl2fch
小程序企图欺骗之,那样即使不报错也可能存在基组数据或轨道不对应的错误。
致谢
感谢wzkchem5,wsr和zhigang的修改建议。 MOKIT中文简介
代码语言:javascript复制https://gitlab.com/jxzou/mokit/-/blob/master/README_zh.md
MOKIT在线手册
代码语言:javascript复制https://jeanwsr.gitlab.io/mokit-doc-mdbook/
引用MOKIT的发表文章一览
代码语言:javascript复制https://jeanwsr.gitlab.io/mokit-doc-mdbook/citing.html