自动做多参考态计算的程序MOKIT

本公众号之前发过几篇多组态（multi-configurational）方法的介绍：

• 用Gaussian做CASSCF计算
• 用Gaussian寻找圆锥交叉点
• 广义价键波函数（GVB）简介
• 广义价键计算及初始轨道的构造

GVB和CASSCF在精度上通常仅是定性正确，定量上差强人意，若需要高精度的结果或与实验值对比，需要进一步做多参考态（multi-reference）方法的计算。GVB、CASCI、CASSCF和DMRG这些（本质上）是多组态波函数，它们是多参考态方法的参考态（reference）。最最常见的多参考态方法有CASPT2、NEVPT2和MRCISD三种，基于CAS或DMRG参考态的都有很多文章发表，相应可以称为CASSCF-NEVPT2、DMRG-NEVPT2等。D. G. Truhlar和L. Gagliardi等人还提出过基于CAS的MC-PDFT方法，后来也推广到了DMRG-PDFT，由于动态相关是用DFT考虑的，比前述几种多参考态方法计算上经济一些。

多组态或多参考态方法计算步骤复杂，要求用户有丰富的计算经验和较强的化学直觉，很多计算（注意仅是写输入文件就）需要有电子结构方法的基础知识。近年来有不少半自动或全自动做多参考态计算的文章发表，意图使这些计算像HF/DFT计算一样简便，但是基本是在文献上或某些课题组里，可获取的程序极少。

笔者前段时间将这几年有关课题中写过的小程序完善和整理了一下，发布了一个自动做多参考态计算的MOKIT程序，顾名思义，分子轨道工具箱。MOKIT含有常见量化软件间传轨道的小程序（见下图）

这些小程序考虑了各个量化程序的基函数的顺序问题（最高到H角动量）、重叠积分问题等等，可以十分精确地在各个量化程序间传轨道，两个程序中同一个波函数方法（例如R(O)HF、UHF、CASSCF等）的电子能量相差小于10-6 a.u.，传轨道过去后通常1圈即收敛。

在此基础上，MOKIT提供了自动做多参考态计算的automr小程序。它可以自动调用这些传轨道的小程序完成系列复杂的操作，如在高斯里算HF，到GAMESS里算GVB，再到下一个程序里算CASSCF，最后到另一个程序里算NEVPT2。MOKIT拥有很多多参考态方法的接口，也列在上图中了。经常有小伙伴辛辛苦苦在一个软件里做了初步计算，换了另一个软件还得从头算，无法利用前一个软件算完的结果，有时候收敛出的能量还不一样。而使用MOKIT则没有这些问题，以下是MOKIT的安装和使用介绍。

1. 下载和安装MOKIT

到https://gitlab.com/jxzou/mokit上下载最新版MOKIT程序。若有读者不想做多参考计算，只想使用传轨道的小程序，可直接下载预编译的Windows版

https://gitlab.com/jxzou/mokit/-/releases

若需要做多参考计算，请下载整个源码压缩包，（发送到Linux系统下再）解压

unzip mokit-master.zip
mv mokit-master mokit

进入文件夹编译

cd mokit/src
make all

耗时约2 min。完成后在~/.bashrc里写上环境变量

export MOKIT_ROOT=/home/$USER/software/mokit
export PATH=$MOKIT_ROOT/bin:$PATH
export PYTHONPATH=$MOKIT_ROOT/lib:$PYTHONPATH
export ORCA=/home/$USER/software/orca-4.2.1/orca
export GMS=/home/$USER/software/games/rungms

以上路径仅为示例，请根据自己的实际情况修改。其中变量ORCA和GMS对应量化软件ORCA和GAMESS可执行文件的完整路径。而PySCF、OpenMolcas和Molpro软件分别由python、pymolcas和molpro命令运行，无需告诉MOKIT它们的位置。除此之外，MOKIT会识别系统中的GAUSS_EXEDIR变量，必要时自动调用Gaussian软件（无论g03/g09/g16），因此也无需知晓高斯具体安装位置。

编译MOKIT需要Fortran编译器（默认ifort）和f2py编译器，运行时还需要一些基本的python库。笔者推荐安装Intel编译器全家桶和Anaconda Python，省事。若想使用gfortran编译器，请自行打开Makefile文件将前几行gfortran相关注释激活（去掉#号），并注释ifort相关语句。

在运行automr前我们还需修改GAMESS源代码。这是因为官方GAMESS只支持少于12对的GVB计算，现今借助于automr可以很容易、自动地做几十对甚至上百对的GVB计算。读者无需手动修改代码，只需将mokit/src/目录下的modify_GMS1.sh和modify_GMS2.f90两个文件拷贝到gamess/source/目录下，运行./modify_GMS1.sh即可。程序包中doc目录下的MOKIT_manual.pdf文件即为程序手册。

前段时间的帖子《广义价键计算及初始轨道的构造》中介绍的自动做GVB计算的代码，其中传轨道部分用的是MOKIT里两年前的一些小程序，且仅支持Cartesian型基函数，有不少限制。如今应当使用MOKIT来做自动化的GVB计算，它默认使用球谐型基函数，减少可能的线性相关问题，自动构造初始轨道的算法也有所升级。

2. 如何引用

若您在您的研究中使用了MOKIT的任何一个子程序或模块，请引用

Molecular Orbital Kit (MOKIT), https://gitlab.com/jxzou/mokit (accessed month day, year)

若您使用MOKIT进行了GVB方法的相关计算，请引用以下2篇文献

DOI: 10.1021/acs.jctc.8b00854; DOI: 10.1021/acs.jpca.0c05216

除此之外，使用者还需引用被automr程序调用到的量化软件。更详细的引用说明请参见程序pdf手册1.2部分。

3. 使用示例

假设使用者已经安装好了常见的量化软件。若未安装，可参考本公众号发过的安装教程：

• Linux下Gaussian 16安装教程
• ORCA软件安装教程
• GAMESS编译教程
• 离线安装PySCF程序（1.5及更高版本）
• OpenMolcas 与 QCMaquis 的安装
• Block-1.5的编译和安装
• Boost.MPI的编译
• 安装基于openmpi的mpi4py

不是所有量化软件都会在一次计算中被调用到的。例如要做NEVPT2计算，automr默认调用PySCF程序做NEVPT2。若用户指定关键词NEVPT2_prog=molpro，则显然需要安装Molpro软件。总的来说，用户应当至少安装Gaussian、PySCF和GAMESS这三个量化软件。automr程序默认调用的都是每个量化软件最快或最稳健的功能，一般无需更改。

程序包中的mokit/examples/目录下提供了许多简单的示例，方便用户练习。注意其中的基组未必合理，有些是为了简化计算量用的。发表级别的计算，应至少采用TZ级别的基组（如cc-pVTZ、def2-TZVP、TZVP等）。以下举几个例子

（1）水分子的CASSCF(4,4)/cc-pVDZ计算，文件名00-h2o_cc-pVDZ_1.5.gjf

%mem=4GB
%nprocshared=4
#p CASSCF/cc-pVDZ

mokit{}

0 1
O   -0.23497692    0.90193619   -0.068688
H    1.26502308    0.90193619   -0.068688
H   -0.73568721    2.31589843   -0.068688

输入文件与Gaussian的格式一样，没有额外的学习成本。automr会自动调用各个量化程序做HF、GVB和CASSCF等计算，过程中会自动确定活性空间为CAS(4,4)，无需用户人为指定。想做GVB、NEVPT2、CASPT2、MRCISD或MCPDFT计算的话，把#p那行中方法名换成对应的方法就行。无需写%chk，因为没有用。注意必须在Title Card那一行写上mokit{}。执行

automr 00-h2o_cc-pVDZ_1.5.gjf>& 00-h2o_cc-pVDZ_1.5.out &

即提交该计算。输出文件末尾展示如下

…
GVB(4)

E(GVB) =       -75.91325105 a.u.
Leave subroutine do_gvb at Tue Dec 29 16:54:33 2020


Enter subroutine do_cas...
DMRG-CASSCF(4,4) using program pyscf
doubly_occ=   3    nvir=  17
No. of active alpha/beta e = 2/2

E(CASCI)  =       -75.90802166 a.u.
E(CASSCF) =       -75.90823032 a.u.
Leave subroutine do_cas at Tue Dec 29 17:00:50 2020

Normal termination of AutoMR at Tue Dec 29 17:00:50 2020

无论调用哪个量化程序做CASSCF计算，automr最终都会给出*CASSCF_NO.fch文件，存有CASSCF自然轨道和轨道占据数，可以用GaussView或Multiwfn+VMD联用查看轨道形状。对自然轨道不了解的小伙伴可以阅读

• 从密度矩阵产生自然轨道-理论篇
• 从密度矩阵产生自然轨道_实战篇（上）

传统做法是局限在某个量化软件中做CASSCF计算，往往只能用其特定的可视化软件看轨道，十分受限。

（2）提供fch(k)文件给automr读入，文件名05-N2_cc-pVQZ_6D10F_4.0.gjf

%mem=4GB
%nprocshared=4
#p NEVPT2/cc-pVQZ

mokit{readuhf='N2_cc-pVQZ_6D10F_4.0_uhf.fchk',ist=1,cart}

fch文件名需写在一对单引号中。由于提供了fch文件，底下无需写电荷、坐标和自旋多重度等信息。mokit{}中的关键词readrhf/readuhf/readno分别用于读取含RHF/UHF/自然轨道的fch文件。有时用户事先通过片段组合波函数等复杂计算得到了一个特定UHF解的文件，想用于多参考态计算，可以用此功能读取。

关键词ist=1表示使用第一种策略，若使用不同的策略（有的）无需GVB计算，详情请阅读pdf手册。关键词cart表示使用Cartesian型基函数（程序默认是球谐型基函数，即5D 7F）。请记得在获得.fch(k)文件前在高斯输入文件里加上关键词nosymm int=nobasistransform。

automr会在计算中自动判定活性空间为CAS(8,8)。十分有经验的用户可以一开始就在输入文件中指定NEVPT2(8,8)，一般无需指定（除非活性空间与预期不同）。另外，若计算中检测到活性空间超过(14,14)，automr会将CASSCF自动切换为DMRGSCF，CASSCF-NEVPT2自动切换为DMRG-NEVPT2（假定用户已经安装好了PySCF和Block程序，否则到这步会报错）。

4. 当前程序的限制

没有万能的程序。MOKIT当前有一些限制，例如

（1）无法直接进行激发态的多参考计算，不含自动产生激发态活性空间的算法。但用户可以基于automr算完产生的文件继续（手动）做激发态计算。

（2）无法使用任何对称性。通常的分子没什么对称性，但对高对称性分子无法节省计算时间。

（3）只能算单点，无法进行结构优化、过渡态、搜索MECP等复杂任务。但是支持在mokit{}中写force关键词输出CASSCF解析梯度。用户也可以用算完的文件（手动）添加关键词做结构优化。

5. 可能出现的问题

当遇到MOKIT使用问题时，请先仔细阅读报错信息，并结合程序手册检查自己输入文件是否写错。手册附录A1部分列出了一些常见问题。另外，应去GitLab下载、安装最新版MOKIT，看问题是否仍存在。

若发现程序bug或想提出建议时，可以向笔者反映（必要时请将输入、输出文件压缩后附上）。有如下途径：

(1) 到https://gitlab.com/jxzou/mokit/-/issues提出问题；

(2) 通过电子邮件njumath[at]sina.cn联系笔者；

(3) 加入MOKIT用户交流QQ群，群号：470745084。

当然，作者不能保证及时回复。。。若遇到量化软件安装问题（非MOKIT问题），请阅读上述3中列出的安装教程，还有问题可以到对应软件的论坛上去提问。若研究的体系十分复杂，或有实验组看到此文，欢迎找笔者进行合作。

6. 注意事项

（1）若提供.fch(k)文件给automr程序读取，请记得在获得.fch(k)文件前在高斯输入文件里加上关键词nosymm int=nobasistransform，避免后续产生不必要的错误。

（2）若使用各个小程序传DFT轨道，由于各个量化程序的积分格点、泛函定义等不尽相同，因此经常无法1圈收敛，但也能在几圈内很快收敛。

（3）尽管MOKIT程序的计算过程是全自动的，无需人为干预。但用户仍需具备使用常见量子化学软件的基本技能（例如熟悉Gaussian软件的常规DFT计算）。若使用者是一名量化新手，应先学习并熟练使用Gaussian软件做常规计算，否则很可能难以正确理解MOKIT的输出内容，或做出错误解读。特地强调这点是因为有些老师看到公众号博文不分难易就转发给小白学生。。。

（4）automr程序输入文件无需人为指定活性空间，但使用者在发表文章时必须写出活性空间大小，及初始轨道如何得到的。不讲明这些细节的多参考态计算是没有意义的。