在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例

在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例

The way of generating natural orbitals based on fch file in Multiwfn and analysis instances about excited state wavefunctions and spin natural orbitals

文/Sobereva @北京科音

First release: 2018-Jan-8  Last update: 2020-Nov-6
 
  

0 前言

在《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》(http://sobereva.com/379)和Multiwfn手册第四章开头都说过,对于Gaussian用户,如果要在Multiwfn中分析后HF波函数、TDDFT激发态波函数或者双杂化泛函的波函数,需要产生相应级别的自然轨道并让Multiwfn载入。自然轨道可以让Gaussian以.wfn/wfx形式输出,也可以让Gaussian将之存到.fch里。前者很容易,但此类输入文件在Multiwfn里没法进行需要利用基函数信息的分析(如Mayer键级、Mulliken电荷);后者虽然没这个限制,但是导出过程很繁琐,需要Gaussian中运行两步,要用的关键词容易忘,而且之后还得在.fch开头自行填上saveNO。

fch文件里是包含密度矩阵信息的,自然轨道也正是通过对角化密度矩阵得到的。但fch里的密度矩阵不会被Multiwfn所直接利用,Multiwfn从输入文件里只会读取轨道信息。从2018-Jan-7更新的Multiwfn开始,为了方便用户,笔者在Multiwfn中添加了一个功能,可以读取fch里的密度矩阵并产生各类自然轨道。本文目的是说一下相关细节,演示一下怎么用,并顺带简要讲一下激发态波函数的分析和自旋自然轨道的分析。

 

1 关于fch里的密度矩阵和自然轨道

在chk里,以及转换出的fch文件里,至少包含了SCF级别的密度矩阵,对于后HF、双杂化泛函、TDDFT等任务,还会包含其它类型密度矩阵(前提是用了density或某些特殊关键词,而且当前理论方法在Gaussian中有解析梯度)。在fch文件里搜索Density就能知道包含哪些。这里罗列一下常见情况:
HF、DFT、CASSCF任务:SCF
CIS、TDDFT任务:CI、CI Rho(1)
CCSD任务:CC
MP2、双杂化泛函任务:MP2
MP3任务:MP3
MP4SDQ任务:MP4
如果你的体系是开壳层(对于CIS/TD来说指的是参考态),还会包含自旋密度矩阵。比如,对一个三重态体系,如果你用的是比如# MP2/cc-pVDZ density关键词,则fch文件里会有这几个字段:
Total SCF Density:HF级别的总密度矩阵
Spin SCF Density:HF级别的自旋密度矩阵
Total MP2 Density:MP2级别的总密度矩阵
Spin MP2 Density:MP2级别的自旋密度矩阵
如果体系是闭壳层,就没有Spin SCF Density和Spin MP2 Density了。

所谓总密度矩阵,就是指Alpha和Beta密度矩阵的加和;而自旋密度矩阵,就是Alpha密度矩阵减去Beta密度矩阵。

Multiwfn能产生的自然轨道有三类,和密度矩阵的关系如下:
(a)自然轨道(Natural orbital, NO):对总密度矩阵对角化得到。本征值是占据数,原则上范围是0.0~2.0,本征矢就是自然轨道系数。这种自然轨道也被称为Spatial natural orbital,因为不含自旋信息。如果体系是开壳层的,这种自然轨道也被确切叫做unrestricted natural orbital (UNO)。
(b)Alpha自然轨道、Beta自然轨道:分别对Alpha密度矩阵和Beta密度矩阵对角化得到,分别描述Alpha和Beta电子,原则上占据数在0.0~1.0之间。这类自然轨道也被称为自然自旋轨道(natural spin orbital, NSO)。
(c)自旋自然轨道(Spin natural orbital, SNO):对自旋密度矩阵对角化得到,占据数在-1.0至1.0之间。占据数为正和为负的SNO分别描述Alpha单电子和Beta单电子。

 

2 实例1:分析甲醛TDDFT计算的S1态波函数

运行以下输入文件,对甲醛在TD-PBE0/6-31G*级别下做电子激发计算,并且把S1态密度矩阵存到chk里(因为TD里的root选项默认为1)。
%chk=C:\H2CO.chk
# PBE1PBE/6-31G* TD density

B3LYP/6-31G* opted

0 1
 C                  0.00000000    0.00000000   -0.56221066
 H                  0.00000000   -0.92444767   -1.10110537
 H                 -0.00000000    0.92444767   -1.10110537
 O                  0.00000000    0.00000000    0.69618930

用formchk将chk转换为fch,里面此时会有Total SCF Density、Total CI Rho(1) Density、Total CI Density三个字段,分别对应于基态DFT密度矩阵、TDDFT非弛豫的S1态密度矩阵(可无视)、TDDFT弛豫的S1态密度矩阵。

下面我们来考察一下S1态的C-O Mayer键级以及O的ADCH原子电荷。启动Multiwfn,依次输入
H2CO.fch
200
16  //基于fch里的密度矩阵产生自然轨道
CI  //载入的密度矩阵级别
由于输入的是CI,所以Multiwfn从fch中读取了Total CI Density字段记录的S1态总密度矩阵,并立刻产生了相应的S1态自然轨道,占据数直接输出在屏幕上了(占据数稍微超过理论上0.0~2.0范围是正常情况,不用纠结):
Occupation numbers:
   2.033734    2.002238    2.000812    2.000459    2.000148    2.000000
   2.000000    0.996739    0.970938    0.000094    0.000062    0.000049
   0.000019    0.000005    0.000002    0.000000    0.000000    0.000000
   0.000000    0.000000   -0.000000   -0.000000   -0.000000   -0.000000
  -0.000000   -0.000000   -0.000000   -0.000000   -0.000238   -0.000366
  -0.000454   -0.000828   -0.001168   -0.002245
此时,内存里的基函数信息已经被更新成了对应S1自然轨道的情况了,做各种基于基函数的分析的结果将是对应S1态的了。我们还要计算ADCH电荷,这是要用GTF信息的,但目前GTF系数还是对应基态DFT轨道的的(因为一开始载入的.fch里的轨道是基态DFT的),且现在轨道占据数已经变成了自然轨道占据数,所以之后做基于GTF信息的分析将没有任何物理意义。为了能计算S1态的ADCH电荷,我们此时应当选y,让程序把当前的波函数信息导出为当前目录下的new.mwfn,并自动载入之,此时内存里的GTF信息也变成对应S1态的了。这个new.mwfn文件我们以后也再可以反复利用。

接着输入以下内容
0  //退回主菜单
9  //键级计算
1  //Mayer键级
得知C-O键级为1.18。然后输入
n  //不导出键级矩阵
0  //退回主菜单
7  //布居分析
11  //ADCH电荷
1  //用内置的球对称化的自由原子密度
得知O的ADCH电荷为-0.128。

我们可以对比一下基态的情况。重启Multiwfn,载入H2CO.fch,此时内存里的基函数信息和GTF信息都是对应基态的,重做上述计算,发现C-O Mayer键级为1.98,O的ADCH电荷为-0.296。

之所以激发态的C-O键级比基态的小得多,激发态的O的电荷也比基态的正得多,是因为此激发是n->pi*激发,一方面削弱了C-O键,一方面一部分电子从氧转移给了C。


3 实例2:分析丁烷双自由基的自旋自然轨道

丁烷双自由基体系在《CASSCF计算双自由基以及双自由基特征的计算》(http://sobereva.com/264)中做了很详细的讨论,建议看看。本节的例子探究一下此体系的SNO轨道。

运行以下输入文件
%chk=C:\C4H8.chk
#p UB3LYP/def2SVP guess=mix nosymm

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


将chk转换成fch后,启动Multiwfn,依次输入
C4H8.fch
200
16
SCF   //考察DFT级别的密度矩阵

此时程序问你,是产生空间自旋轨道(此时叫做UNO),Alpha和Beta各自的自然自旋轨道(NSO),还是自旋自然轨道(SNO)。我们选3产生SNO,看到的占据数如下
    0.957908    0.045983    0.043627    0.037531    0.021487    0.021280
    0.020200    0.017324    0.008494    0.007915    0.006743    0.006539
    0.000626    0.000604    0.000082    0.000082    0.000000    0.000000
...略
   -0.000000   -0.000000   -0.000082   -0.000082   -0.000604   -0.000626
   -0.006539   -0.006743   -0.007915   -0.008494   -0.017324   -0.020200
   -0.021280   -0.021487   -0.037531   -0.043627   -0.045983   -0.957908
可见占据数为正的SNO当中只有一个数值较大的,即第一个(0.9579);而占据数为负的SNO当中也只有一个数值较大的,即最后一个(-0.9579)。而且,我们也知道,由于此体系是个很理想的双自由基体系,本来就应该有一个Alpha和一个Beta单电子。因此,未成对的Alpha电子和Beta电子分布特征通过考察这两个轨道就足够得知了。

我们选y产生new.mwfn并自动载入之,然后进入主功能0观看这两个轨道的特征。值得一提的是,new.mwfn被Multiwfn视为是开壳层的,因此形式上会有Alpha和Beta两套轨道,而我们产生的SNO轨道只有一套,此时SNO轨道占用的是用于记录Alpha轨道那部分空间(绝对不是说SNO的物理意义是Alpha轨道),而Beta轨道那部分没有意义,占据数也都是0。我们在主功能0的界面里可以选左上角的Orbital info.,从文本窗口里找出占据数为0.957908和-0.957908的轨道,可知序号分别是1和96(其实96就是体系的基函数数目,因此都没必要去列表里找,直接选就行),这两个轨道的图形如下:



可见,Alpha单电子主要分布在C1上,在C4-C7之间也有少量分布。Beta单电子主要分布在C10上,同样在C4-C7之间也有少量分布。

如果我们想定量知道比如Alpha单电子在各个原子上的量,我们可以对主要表现Alpha单电子的SNO1做轨道成分分析。我们回到主菜单,依次输入
8  //轨道成分分析
8  //Hirshfeld方式做轨道成份分析
1  //用内置的球对称化的自由原子密度
1  //1号轨道
结果为
Atom     1(C ) :   70.148760%
Atom     2(H ) :    6.317230%
Atom     3(H ) :    6.317230%
Atom     4(C ) :    7.588237%
Atom     5(H ) :    1.254212%
Atom     6(H ) :    1.254212%
Atom     7(C ) :    5.065489%
Atom     8(H ) :    0.591654%
Atom     9(H ) :    0.591654%
Atom    10(C ) :    0.758644%
Atom    11(H ) :    0.056339%
Atom    12(H ) :    0.056339%
数值和我们从轨道图形上看到的很好相符。

众所周知,研究单电子分布最常用的方法是分析自旋密度。我们重新载入C4H8.fch,绘制一下自旋密度图,步骤见《谈谈自旋密度、自旋布居以及在Multiwfn中的绘制和计算》(http://sobereva.com/353),图像如下:


我们再用主功能15,在Hirshfeld划分下计算自旋密度在各个原子空间中的积分值(原子自旋布居):
Atomic space        Value
  1(C )            0.70928760
  2(H )            0.04488357
  3(H )            0.04488357
  4(C )           -0.01371421
  5(H )            0.00849956
  6(H )            0.00849956
  7(C )            0.01371421
  8(H )           -0.00849956
  9(H )           -0.00849956
 10(C )           -0.70928760
 11(H )           -0.04488357
 12(H )           -0.04488357

无论是从等值面图上,还是原子自旋布居上,都看到自旋密度分布特征基本等于主要描述Alpha单电子的SNO1和主要描述Beta单电子的SNO96特征的叠加。对于两个SNO基本没有重叠的体系两端的亚甲基的碳,其自旋布居和它对SNO的贡献量十分接近,而在C4、C7部分,由于两个SNO重叠厉害,导致显著的相互抵消,所以它们虽然对SNO贡献不很小,但是自旋布居却非常小。

可见,分析SNO的好处是可以分别考察未成对的Alpha和Beta电子原本是怎么分布的,且可以从轨道分析的角度研究;而考察自旋密度,则便于了解在Alpha和Beta单电子分布发生一定程度抵消后的单电子“净”分布情况,分析也更为省事,不用讨论多个轨道。两种分析是互补的。