使用Multiwfn计算激发态间的跃迁偶极矩和各个激发态的偶极矩

使用Multiwfn计算激发态间的跃迁偶极矩和各个激发态的偶极矩

文/Sobereva @北京科音
First release: 2014-Mar-28   Last update: 2022-Apr-13


1 前言

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

强大的波函数分析程序Multiwfn的电子激发分析模块中的子功能5能够计算所有态(包括基态和激发态)之间的跃迁电和磁偶极矩,可以用于上述研究。此外,这个功能还可以直接给出各个激发态的电偶极矩,通过与基态的偶极矩求差就可以考察电子激发过程中偶极矩的变化。在本文就简要演示一下这个功能的使用。请读者务必使用官网上最新版Multiwfn,可以在http://sobereva.com/multiwfn免费下载。如果对Multiwfn一无所知的话,建议先阅读《Multiwfn入门tips》(http://sobereva.com/167)和《Multiwfn FAQ》(http://sobereva.com/452)。

对于Gaussian用户,使用Multiwfn的上述功能需要两类文件:(1)Gaussian的CIS、TDHF、TDDFT或TDA-DFT的输出文件 (2)相应任务产生的.fch/fchk文件。Multiwfn的这个功能不限于Gaussian用户使用,也支持ORCA、GAMESS-US等程序,详见Multiwfn手册3.21.A节的说明。


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 实例

下面例子涉及的文件都可以在http://sobereva.com/attach/227/file.rar下载。

下面以苯酚为例进行说明怎么利用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输出到当前目录下的transdipmom.txt里。程序先给出了基态的偶极矩,然后给出基态到各个激发态的跃迁偶极矩、激发能和振子强度:
 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 the case of i=j corresponds to contribution of electron t
o dipole moment of excited state i
 Transition electric dipole moment between excited states (a.u.):
     i     j         X             Y             Z        Diff.(eV)   Oscil.str
     1     1    -0.5515307     0.6638909     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.5370763     0.7195150     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.0119437    -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.2932424     2.0434979     0.0000000     0.00000     0.00000
     4     5     0.0000000     0.0000000     0.0048938     0.22550     0.00000
     5     5    -0.4634208     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给出的2.0119437    -3.9241936     0.0000000很相符。注意这里用了rhoci关键词,如果只写density的话,那么Gaussian传递给L601模块的激发态的密度是弛豫的密度。而当用rhoci时,传递的是非弛豫的密度,这才是和通过组态系数和基态轨道直接产生的密度直接对应的,也正是与Multiwfn输出的是对应的。

如果你想得到各个态之间的跃迁磁偶极矩,进入功能18的子功能5后先选0,选Magnetic,然后再按上面做法操作即可。注意跃迁磁偶极矩定义为-i<i|r×▽|j>,而Multiwfn输出的值对应的是<i|r×▽|j>部分,这和Gaussian输出文件里给出的跃迁磁偶极矩用的形式相同。

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

 Excited state electric dipole moments (a.u.):
  State         X             Y             Z        exc.(eV)    exc.(nm)
     1     -0.551525      0.663891      0.000000      5.1144      242.42
     2     -0.537071      0.719515      0.000000      5.9365      208.85
     3      2.011949     -3.924194      0.000000      6.1432      201.82
     4     -0.293237      2.043498      0.000000      6.6769      185.69
     5     -0.463415      1.058643      0.000000      6.9024      179.62


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