将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法

将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法

文/Sobereva@北京科音  2019-Oct-11


1 说明

ORCA是非常强大的量子化学程序,笔者之前也写过不少相关文章,开发了不少辅助工具,见sobereva.com右侧的ORCA分类。ORCA相对于最常用的量子化学程序Gaussian来说,有一个缺点是SCF收敛做得不够好,很多Gaussian能收敛的情况ORCA难以收敛,而且Gaussian的SCF不收敛的解决方案比较成熟,见《解决SCF不收敛问题的方法》(http://sobereva.com/61)。另外,ORCA在产生初猜波函数方面也没Gaussian灵活,比如Gaussian用户可以用GaussView构建片段组合波函数作为初猜,而且利用stable=opt关键词还可以自动优化出稳定的波函数,这对于双自由基、反铁磁性耦合体系的研究十分重要,参看《谈谈片段组合波函数与自旋极化单重态》(http://sobereva.com/82)。虽说ORCA也不是没法算这类自旋极化单重态体系,利用%SCF FlipSpin也可以实现,但往往明显不如Gaussian方便。

显然,如果能让其它量子化学程序,特别是波函数初猜以及SCF迭代方面做得很好的Gaussian产生的波函数作为ORCA的初猜波函数,就能令ORCA取长补短,从很多恼人的问题中解脱。

Multiwfn(主页&下载地址:http://sobereva.com/multiwfn)程序支持载入.fch、.molden、GAMESS-US/Firefly输出文件这些含有基组定义以及轨道信息的格式(在Multiwfn里称为“基函数信息”),并且可以导出各种含有波函数信息的格式,支持的文件格式详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》(http://sobereva.com/379)。Multiwfn支持绝大多数主流量化程序产生的波函数,并可以作为波函数文件格式的转换器来用。从2019-Oct-10更新的Multiwfn开始,Multiwfn的导出文件功能更是新增了产生.mkl文件的功能。.mkl文件是老版本Molekel的输入文件,可以被ORCA自带的orca_2mkl工具转换成.gbw文件,ORCA计算时可以从.gbw文件中读取轨道作为初猜波函数。显然,利用Multiwfn可以直接实现将其它量化程序产生的波函数作为ORCA的初猜波函数的目的。甚至可以说,只要存在一种量化程序有办法得到当前体系当前级别下收敛的波函数,用ORCA计算也一定能收敛。

有一点要提及的是ORCA是基于球谐型高斯函数做的计算,因此用Multiwfn实现这个目的,其它程序计算时也应当用球谐型高斯函数。不了解这方面的话看《谈谈5d、6d型d壳层基函数与它们在Gaussian中的标识》(http://sobereva.com/51)。


2 例子:丁烷双自由基

这里举个例子,计算丁烷双自由基C4H8,用以说明如何将Gaussian收敛的波函数作为ORCA的初猜波函数,请大家举一反三。本例涉及的文件都在http://sobereva.com/attach/517/file.rar里。

本例使用2019-Oct-10更新的Multiwfn,Gaussian使用G16 A.03版,ORCA使用4.2版。

首先用Gaussian运行C4H8.gjf,内容如下。任务是对C4H8在B3LYP/def2-SVP级别下产生最稳定波函数,对此体系也是对称破缺波函数。使用def2-SVP的时候Gaussian默认用的是球谐型基函数。

%chk=C:\C4H8.chk
# B3LYP/def2SVP stable=opt
[空行]
ub3lyp/6-31g(d) opted
[空行]
0 1
 C                 -0.74400100    1.78566400    0.00000000
 H                 -0.60282700    2.33865300    0.92499500
 H                 -0.60282700    2.33865300   -0.92499500
 C                 -0.74400100    0.30988100    0.00000000
 H                 -1.25452600   -0.08746700    0.88463900
 H                 -1.25452600   -0.08746700   -0.88463900
 C                  0.74400100   -0.30988100    0.00000000
 H                  1.25452600    0.08746700   -0.88463900
 H                  1.25452600    0.08746700    0.88463900
 C                  0.74400100   -1.78566400    0.00000000
 H                  0.60282700   -2.33865300   -0.92499500
 H                  0.60282700   -2.33865300    0.92499500

然后用formchk将chk转换为fch,载入Multiwfn后依次输入
100  //其它功能 Part 1
2   //导出文件
9   //导出mkl文件
C4H8.mkl
现在当前目录下就有了C4H8.mkl。

在当前目录下运行orca_2mkl C4H8 -gbw,就把C4H8.mkl转换为了C4H8.gbw。下面,我们将C4H8.gbw里的轨道作为初猜波函数,用ORCA也在B3LYP/def2-SVP下跑一下这个双自由基的单点。输入文件如下,将之命名为C4H8.inp,把它和C4H8.gbw都放在当前目录下,然后进行计算。运行时ORCA会自动从C4H8.gbw中读取轨道作为初猜波函数。由于C4H8.gbw里的波函数是非限制性波函数,因此写了UKS关键词。由于ORCA的B3LYP和Gaussian的定义不同,所以加了/G来和Gaussian一致。此文件里的坐标和上面Gaussian输入文件里的坐标精确一致。

! UKS B3LYP/G def2-SVP miniprint nopop
* xyz   0   1
 C                 -0.74400100    1.78566400    0.00000000
 H                 -0.60282700    2.33865300    0.92499500
 H                 -0.60282700    2.33865300   -0.92499500
 C                 -0.74400100    0.30988100    0.00000000
 H                 -1.25452600   -0.08746700    0.88463900
 H                 -1.25452600   -0.08746700   -0.88463900
 C                  0.74400100   -0.30988100    0.00000000
 H                  1.25452600    0.08746700   -0.88463900
 H                  1.25452600    0.08746700    0.88463900
 C                  0.74400100   -1.78566400    0.00000000
 H                  0.60282700   -2.33865300   -0.92499500
 H                  0.60282700   -2.33865300    0.92499500
 *

迭代过程信息如下
ITER       Energy         Delta-E        Max-DP      RMS-DP      [F,P]     Damp
               ***  Starting incremental Fock matrix formation  ***
  0   -157.0020361024   0.000000000000 0.00099750  0.00002892  0.0004319 0.7000
  1   -157.0020389471  -0.000002844760 0.00101224  0.00002810  0.0003439 0.7000
                               ***Turning on DIIS***
  2   -157.0020411727  -0.000002225540 0.00299633  0.00007975  0.0002632 0.0000
  3   -157.0020467355  -0.000005562801 0.00029108  0.00000721  0.0000572 0.0000
  4   -157.0020467759  -0.000000040442 0.00006147  0.00000180  0.0000542 0.0000
                 **** Energy Check signals convergence ****

               *****************************************************
               *                     SUCCESS                       *
               *           SCF CONVERGED AFTER   5 CYCLES          *
               *****************************************************

可见由于用了Gaussian收敛的波函数当做初猜,收敛非常快和顺利,第一轮迭代时的能量和密度矩阵变化就已经极小了,之后很快就收敛了。之所以不是一轮就收敛,是因为Gaussian和ORCA用的DFT积分格点有异。如果二者都用的是HF方法的话,ORCA里仅仅一轮就能收敛。

如果gbw和输入文件不同名,为了能从gbw中读取初猜波函数,应当写上moread关键词,然后加上一行诸如%moinp "/sob/Azusa_Nakano.gbw"来指明读取的gbw文件位置。

顺带一提,利用Multiwfn还能将其它量化程序产生的波函数作为GAMESS-US的初猜波函数,因为如果载入的波函数文件含有基函数信息,Multiwfn产生的GAMESS-US输入文件里可以直接带初猜信息。Multiwfn也可以将其它量化程序产生的波函数当Gaussian的初猜,因为Multiwfn可以导出fch格式,用Gaussian的unfchk转成chk后,就可以用guess=read来从中读取波函数当初猜了。若有不清楚的,参看前述的《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》。