使用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是笔者的个人喜好，但做实际计算至少应采用TightSCF。defgrid3表示使用尽可能高的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

注意事项

目前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内容如下

%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

后续可以使用fch2inp，fch2inporb，fch2com，bas_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