使用MOKIT将ORCA的轨道信息传输到其他量化程序

发表时间: 2023-08-21 09:21

之前我们介绍过从Gaussian或PySCF向其他量化程序传轨道

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

本文介绍如何使用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内容如下

%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是笔者的个人喜好,但做实际计算至少应采用TightSCFdefgrid3表示使用尽可能高的COSX格点。该例使用了RIJCOSX加速,但这对轨道系数影响非常小,更何况用了比较精细的格点和严格的收敛限,预期得到的轨道与一个传统的RHF/def2-TZVP计算的轨道极其接近。用ORCA算完后运行orca_2mkl h2o -mkl转化生成h2o.mkl文件,下面我们以该文件为例展示如何传轨道,读者若使用本文功能,请务必顺带阅读文末注意事项

1. ORCA传轨道给AMESP

mkl2amo h2o.mkl

会产生h2o.amo和h2o.aip两个文件,分别含基组和轨道数据,坐标及关键词。为了让AMESP读入轨道,需要运行

a2m h2o.amo

即将h2o.amo文本文件转化为二进制文件h2o.mo,其中a2m是AMESP自带的小程序

2. ORCA传轨道给BDF

mkl2bdf h2o.mkl

会产生H2O.BAS,h2o.inp和h2o.scforb三个文件,分别含基组数据,坐标及关键词,轨道数据。

3. ORCA传轨道给CFOUR

mkl2cfour h2o.mkl

会产生ZMAT, GENBAS和OLDMOS三个文件,注意CFOUR所用文件名是固定的。mkl2cfour目前有个缺点:在书写ORCA输入文件时需要使用CFOUR输出文件里体系的直角坐标,产生的mkl文件才能给mkl2cfour使用。这是因为CFOUR会旋转分子朝向,只有使用CFOUR输出文件里体系的直角坐标,其已经被旋转至程序标准朝向,才往往不会再被旋转了。如果遇到CFOUR的SCF不收敛,显然此时有输出文件,这点倒不成问题。

4. ORCA传轨道给Dalton

mkl2dal h2o.mkl

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

5. ORCA传轨道给GAMESS

mkl2inp h2o.mkl

会产生h2o.inp文件,内含坐标,基组数据和一些简单的关键词。

6. ORCA传轨道给Molpro

mkl2com h2o.mkl

会产生h2o.a和h2o.com文件,前者含Alpha轨道,后者含坐标,基组和关键词。如果是UHF计算,还会产生h2o.b文件(含Beta轨道)。

7. ORCA传轨道给OpenMolcas

mkl2inporb h2o.mkl

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

8. ORCA传轨道给PSI4

mkl2psi h2o.mkl

会产生h2o.A和h2o.inp文件。前者含Alpha轨道,后者含坐标,基组和关键词。如果是UHF类型计算,还会有h2o.B文件(内含Beta轨道)。

9. ORCA传轨道给PySCF

mkl2py h2o.mkl

会产生h2o.fch和h2o.py文件。在运行.py文件时会从.fch文件里读取轨道。

10. ORCA传轨道给Q-Chem

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小程序

mkl2gjf h2o.mkl

会产生h2o.gjf文件,内含坐标、基组信息和轨道信息。注意Gaussian是支持从gjf文件中读取轨道的(关键词guess=cards)。

11.2 使用mkl2fch小程序

11.2.1 适用于常见HF/DFT/CASSCF轨道

mkl2fch h2o.mkl

若检测到当前目录下无h2o.fch文件,则会从零创建h2o.fch;若已存在h2o.fch文件,则仅将轨道信息写入h2o.fch。

11.2.2 适用于UNO、MP2、CCSD、CASSCF等自然轨道

此处以RI-MP2为例,输入文件计算级别部分为

! 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(自旋)自然轨道

这类计算输入文件较复杂,往往涉及到对称破缺单重态和检验波函数稳定性等问题,计算实例见笔者在计算化学公社论坛上的回帖

http://bbs.keinsci.com/thread-29762-1-1.html

这里我们假设已经获得了.mdci.nat文件,接着执行

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

注意事项

  1. 目前ORCA中开启点群对称性做计算只能旋转分子朝向、识别点群和不可约表示,并不能利用对称性加速计算,所以很少使用。是否开启点群对称性做计算不影响本文功能。但传轨道后在目标程序中应关闭对称性(生成的输入文件里写了这些关键词),笔者没有让各个小程序支持点群对称性的意向。

  2. 以上各个小程序都会产生目标程序输入文件,内含坐标和基组数据,不需要用户再去手动书写基组名称。强烈推荐用户使用该文件进行计算,既能免去手写基组的麻烦,也能保证传轨道时对应性更好。

  3. ORCA仅支持球谐型基函数,不支持笛卡尔型基函数,因此传轨道后产生的文件采用的也是球谐型基函数进行计算。

  4. 若ORCA的SCF有上千a.u.的剧烈振荡,很可能是出现了基函数线性依赖导致的,此时即使侥幸收敛了能量也未必靠谱,需要在输入文件里加上

%scf 
sthresh 1e-6
end

此阈值是Gaussian和GAMESS处理基函数线性依赖的默认值,比较靠谱。

  1. 对于ORCA产生的.mp2nat,.mdci.nat和.nbo这些实际上是gbw格式的文件,转化成mkl文件前除了修改文件名外,还有一种办法
orca_2mkl h2o.mp2nat h2o_mp2no.mkl -mkl -anyorbs

就是两个文件名参数都带上扩展名,然后在最后加-anyorbs参数,可省去改名麻烦。

  1. mkl文件有个缺点是不含赝势信息,笔者曾在ORCA官方论坛上发帖建议未来ORCA在mkl文件中增加赝势信息,但ORCA开发者表示目前没兴趣。若读者在计算中使用全电子基组,自然无此问题;若用了赝势,按上文操作产生其他量化程序的文件不会含赝势信息,即使轨道系数正确,SCF也会剧烈振荡。这里笔者推荐一种解决办法:

Step 1. 先用Gaussian算个几秒的任务 这里以二甲基锌(ZnMe2)为例,假设我们对Zn用SDD,gjf内容如下

%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相关文件 执行

fch2mkl ZnMe2.chk

随即产生ZnMe2.fch,ZnMe2_o.inp和ZnMe2_o.mkl文件,现在我们就有了一个含赝势信息的fch文件。ZnMe2_o.mkl含有来自Gaussian的轨道,除非有特殊计算目的,否则一般推荐使用该mkl文件,即运行

orca_2mkl ZnMe2_o -gbw

将mkl转化为gbw文件。我们用ZnMe2_o.inp文件在ORCA中做目标计算,最后把轨道传回已存在的fch文件

orca_2mkl ZnMe2_o -mkl # 假设ZnMe2_o.gbw已被更新,重新产生mkl
cp ZnMe2.fch ZnMe2_o.fch # 备份一下,以防覆盖
mkl2fch ZnMe2_o.mkl # 将轨道传回ZnMe2_o.fch

后续可以使用fch2inpfch2inporbfch2combas_fch2py等小程序传给其他量化程序做计算,且均含赝势信息。mkl2fch检测到有fch文件存在时会实直接使用该文件,而非从零生成,因此避免了赝势缺失问题。注意不能跳过Step 1随便塞一个fch文件给mkl2fch小程序企图欺骗之,那样即使不报错也可能存在基组数据或轨道不对应的错误。

致谢

感谢wzkchem5,wsr和zhigang的修改建议。

https://gitlab.com/jxzou/mokit/-/blob/master/README_zh.md

MOKIT在线手册

https://jeanwsr.gitlab.io/mokit-doc-mdbook/

引用MOKIT的发表文章一览

https://jeanwsr.gitlab.io/mokit-doc-mdbook/citing.html