利用Multiwfn计算Gaussian输出的激发态之间的跃迁偶极矩

利用Multiwfn计算Gaussian输出的激发态之间的跃迁偶极矩

文/Sobereva @北京科音
First release: 2014-Mar-28   Last update: 2019-Jan-20


1 前言

计算激发态之间的跃迁偶极矩有一些用处,比如用于sum-over-states (SOS)方式计算(超)极化率(详见http://sobereva.com/232)、计算激发态到激发态的吸收光谱(用瞬态吸收光谱技术可以测定)。Gaussian的CIS、TDHF、TDA-DFT、TDDFT可以直接给出基态到激发态的跃迁偶极矩。CIS、TDA也提供了alltransitiondensities关键词可以把激发态之间的跃迁偶极矩输出出来,但是对于TD任务则没有这个关键词。虽然Gaussian有个density=transition=(n,m)关键词可以产生n到m激发态的跃迁密度,然而这个关键词在目前的版本中貌似功能不正常,而且光是有跃迁密度也没用,还得利用偶极矩积分才能得到跃迁偶极矩,但Gaussian也并不给出这样的信息。

笔者多次被问及怎么基于Gaussian得到激发态之间的跃迁偶极矩,对于TD任务确实没有可行的办法。于是笔者在Multiwfn多功能波函数分析程序的电子激发分析模块中加入了这一功能,这里就介绍下怎么使用。另外,这个功能还可以直接给出各个激发态的偶极矩。本文用的是Multiwfn 3.6版。如果对Multiwfn一无所知的话,建议先阅读《Multiwfn入门tips》(http://sobereva.com/167

Multiwfn计算激发态间的跃迁偶极矩需要两个文件:(1)Gaussian的CIS、TDHF、TDDFT或TDA-DFT的输出文件 (2)相应任务的.fch/fchk文件。另外,Multiwfn也支持载入ORCA等其它程序的输出文件来产生各个态之间的跃迁偶极矩,见Multiwfn手册3.21.5节

Gaussian的输入文件对于CIS和TDDFT,分别写成类似这样
# B3LYP/6-31+G(d) TD(nstates=10) IOp(9/40=4)
# CIS(nstates=10)/6-31+G(d) IOp(9/40=4)
这里假设算10个激发态。Multiwfn计算时需要利用Gaussian输出的激发态的组态系数,默认情况下只有绝对值大于0.1的系数会被输出出来,而较小的都不输出,这样的话Multiwfn算出的结果将会不太准确。IOp(9/40=x)的含义是将系数大于10^-x的组态都输出出来,因此IOp(9/40=4)会把系数绝对值大于0.0001的组态都输出。x设的越大,输出的组态越多,结果越准确,但是x太大的话Multiwfn的耗时也会非常大,通常x=3或4就够了,精度足够,计算速度也比较快。

在计算之后得到了.out/.log文件,还同时得到了.chk文件。用formchk将.chk转换为.fch/fchk,这个文件里记录了基函数的定义以及基态轨道信息,这是Multiwfn要利用的。


2 实例

下面以苯酚为例进行说明怎么利用Multiwfn计算激发态之间的跃迁偶极矩,输入文件如下,计算后得到phenol.out以及phenol.fch。
%chk=C:\phenol.chk
# b3lyp/6-311G* TD(nstates=5) IOp(9/40=4)

b3lyp/6-311G** opted

0 1
 C                  0.01810200   -1.86802400    0.00000000
 C                  1.23175000   -1.16732400    0.00000000
 C                  1.23175000    0.23407600    0.00000000
 C                  0.01810200    0.93477600    0.00000000
 C                 -1.19554600    0.23407600    0.00000000
 C                 -1.19554600   -1.16732400    0.00000000
 H                  0.01810200   -2.93802400    0.00000000
 H                  2.15839700   -1.70232400    0.00000000
 H                  2.15839700    0.76907600    0.00000000
 H                 -2.12219300    0.76907600    0.00000000
 H                 -2.12219300   -1.70232400    0.00000000
 O                  0.01810200    2.36477600    0.00000000
 H                 -0.88699500    2.68477600    0.00000000

启动Multiwfn,依次输入以下命令,//后面的是注释。
C:\phenol.fch    //先载入电子激发任务产生的fch文件
18   //电子激发分析功能,包含多个子功能。这些功能无与伦比的强大,建议参看《Multiwfn支持的电子激发分析方法一览》(http://sobereva.com/437
5    //计算激发态间的跃迁偶极矩
C:\phenol.out    //电子激发任务的Gaussian输出文件
此时屏幕上首先输出了激发态的汇总信息
 Exc.state#     Exc.energy(eV)     Multi.   MO pairs    Normalization
       1           5.11440           1         1478        0.500002
       2           5.93650           1         1697        0.500000
       3           6.14320           1          919        0.499999
       4           6.67690           1          828        0.499994
       5           6.90240           1         1593        0.499994
其中N_pairs代表这个激发态在输出文件中通过多少组态来表示,IOp(9/40=x)的x越大显然N_pairs也就越大。Multi.是激发态的自旋多重度。Sum coeff.^2越接近理想值说明结果精度越高,对于闭壳层和开壳层情况理想值分别是0.5和1.0。如果偏离理想值比较大,则应该加大x来保证结果精度。

然后选1,激发态间的跃迁偶极矩就输出到了屏幕上。也可以选2输出到文本文件里。程序先给出了基态的偶极矩,然后给出基态到各个激发态的跃迁偶极矩、激发能和振子强度:
 Ground state dipole moment in X,Y,Z:   -0.554241   -0.157977    0.000000 a.u.

 Transition dipole moment between ground state (0) and excited states (a.u.)
     i     j         X             Y             Z        Diff.(eV)   Oscil.str
     0     1    -0.4634428    -0.0387809     0.0000000     5.11440     0.02710
     0     2    -0.0410220     0.4644091     0.0000000     5.93650     0.03161
     0     3     0.0000000     0.0000000     0.0024911     6.14320     0.00000
     0     4     0.0000000     0.0000000    -0.0776262     6.67690     0.00099
     0     5    -1.3520332    -0.2422995     0.0000000     6.90240     0.31905

接下来,输出的是激发态之间的跃迁偶极矩<i|-r|j>、能量差和振子强度。对于i=j的情况来说,其数值<i|-r|i>对应的是第i激发态的偶极矩的由电子贡献的部分(注意这不等于这个态的偶极矩,因为还有另一部分,即原子核电荷对偶极矩的贡献)。
 Note: In below output and for the cases of i=j, only electronic contribution to
 dipole moment of the state is taken into account
 Transition dipole moment between excited states (a.u.):
     i     j         X             Y             Z        Diff.(eV)   Oscil.str
     1     1    -0.5515250     0.6638910     0.0000000     0.00000     0.00000
     1     2     0.1289203     0.0072212     0.0000000     0.82210     0.00034
     1     3     0.0000000     0.0000000    -0.0575947     1.02880     0.00008
     1     4     0.0000000     0.0000000     0.0104765     1.56250     0.00000
     1     5     0.0412045     1.1173067     0.0000000     1.78800     0.05476
     2     2    -0.5370707     0.7195151     0.0000000     0.00000     0.00000
     2     3     0.0000000     0.0000000    -0.0271832     0.20670     0.00000
     2     4     0.0000000     0.0000000     0.0723899     0.74040     0.00010
     2     5     0.7076947    -0.0614123     0.0000000     0.96590     0.01194
     3     3     2.0119493    -3.9241936     0.0000000     0.00000     0.00000
     3     4     0.3712230     1.1487737     0.0000000     0.53370     0.01906
     3     5     0.0000000     0.0000000    -0.0330295     0.75920     0.00002
     4     4    -0.2932368     2.0434980     0.0000000     0.00000     0.00000
     4     5     0.0000000     0.0000000     0.0048938     0.22550     0.00000
     5     5    -0.4634152     1.0586426     0.0000000     0.00000     0.00000

值得一提的是,如果在Gaussian中使用下面这样的关键词,在输出文件末尾就会由L601模块输出第3个激发态的偶极矩
# b3lyp/6-311G* TD(root=3,nstates=5) density=rhoci
给出的结果为(debye)
X=              5.1144    Y=             -9.9743    Z=              0.0000
换算成a.u.单位后,结果为X=2.0119 Y=-3.9238 Z=0.0000,可见和Multiwfn给出的很相符。注意这里用了rhoci关键词,如果只写density的话,那么Gaussian传递给L601模块的激发态的密度是弛豫的密度。而当用rhoci时,传递的是非弛豫的密度,这才是和通过组态系数和基态轨道直接产生的密度直接对应的,也正是与Multiwfn输出的是对应的。

Multiwfn也可以直接给出各个激发态的包含了原子核贡献的偶极矩。也就是进入功能18的子功能5后选择4,结果就会输出到当前目录下的dipmom.txt中,内容为:
 Note: The dipole moments shown below include both nuclear charge and electronic contributions
 Ground state dipole moment in X,Y,Z:   -0.554241   -0.157977    0.000000 a.u.

 Excited state dipole moments (a.u.):
   State      X           Y           Z        exc.(eV)    exc.(nm)
     1   -0.551519    0.663891    0.000000      5.1144      242.59
     2   -0.537065    0.719515    0.000000      5.9365      209.00
     3    2.011955   -3.924194    0.000000      6.1432      201.96
     4   -0.293231    2.043498    0.000000      6.6769      185.82
     5   -0.463410    1.058643    0.000000      6.9024      179.75

PS:之所以这里把原子核对电子态的贡献也考虑后,偶极矩数值和之前只考虑电子贡献的时候一样,那是因为Gaussian默认会把体系原点放置到核电荷中心位置,此时原子核对偶极矩贡献恰好为0。