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

后记:2017-Aug-19及之后更新的Multiwfn也支持载入ORCA输出文件来产生各个态之间的跃迁偶极矩,见Multiwfn手册3.21.5节。


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

文/Sobereva(3)
First release: 2014-Mar-28   Last update: 2016-Aug-31


Gaussian的CIS、TDHF、TDA-DFT、TDDFT可以直接给出基态到激发态的跃迁偶极矩。CIS、TDA也提供了alltransitiondensities关键词可以把激发态之间的跃迁偶极矩输出出来,但是对于TD任务则没有这个关键词。虽然Gaussian有个density=transition=(n,m)关键词可以产生n到m激发态的跃迁密度,然而这个关键词在目前的版本中貌似功能不正常,而且光是有跃迁密度也没用,还得利用偶极矩积分才能得到跃迁偶极矩,但Gaussian也并不给出这样的信息。

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

Multiwfn计算激发态间的跃迁偶极矩需要两个文件:(1)Gaussian的CIS或TDHF或TDDFT的输出文件 (2)相应任务的.fch文件。

Gaussian的输入文件对于CIS、TDHF、TDDFT,分别写成类似这样
# B3LYP/6-31+G(d) TD(nstates=10) IOp(9/40=4)
# CIS(nstates=10)/6-31+G(d) IOp(9/40=4)
# HF/6-31+G(d) TD(nstates=10) 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,这个文件里记录了基函数的定义以及基态轨道信息,这是Multiwfn要利用的。


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

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:\gtest\phenol.fch    //先载入电子激发任务产生的fch文件
18   //电子激发分析功能,包含多个子功能。这些功能无与伦比的强大,建议参看Multiwfn手册4.18节的例子以及3.21节的原理和细节的介绍来学习使用
5    //计算激发态间的所有跃迁偶极矩
c:\gtest\phenol.out    //电子激发任务的Gaussian输出文件
此时屏幕上首先输出了激发态的汇总信息
Exc.state#     Exc.energy(eV)     Multi.   N_pairs    Sum coeff.^2
      1           5.06970           1        449        0.500004
      2           5.47020           1        102        0.499997
      3           5.94190           1         82        0.499975
      4           6.26640           1         80        0.499979
      5           6.38330           1         83        0.499978
其中N_pairs代表这个激发态在输出文件中通过多少组态来表示,IOp(9/40=x)的x越大显然N_pairs也就越大。Multi.是激发态的自旋多重度。Sum coeff.^2越接近理想值说明结果精度越高,对于闭壳层和开壳层情况理想值分别是0.5和1.0。如果偏离理想值比较大,则应该加大x来保证结果精度。

然后选1,结果就输出到了屏幕上。也可以选2输出到文本文件里。结果如下:
Transition dipole moment between excited states (a.u.):
    i     j        X           Y           Z      Ene.diff.(eV)   Oscil.str
    1     1   -0.510929    0.541413    0.000000     0.00000        0.00000
    1     2    0.000000    0.000000    0.180185     0.40050        0.00032
    1     3    0.000000    0.000000   -0.052359     0.87220        0.00006
    1     4    0.000000    0.000000    0.033679     1.19670        0.00003
    1     5    0.000000    0.000000   -0.035481     1.31360        0.00004
    2     2    2.353643   -3.851667    0.000000     0.00000        0.00000
    2     3    0.110576    1.848714    0.000000     0.47170        0.03964
    2     4    0.029210   -0.159245    0.000000     0.79620        0.00051
    2     5    1.825329    0.672088    0.000000     0.91310        0.08464
    3     3   -0.434476    2.840406    0.000000     0.00000        0.00000
    3     4   -0.437467    0.064725    0.000000     0.32450        0.00155
    3     5    3.994461   -0.441695    0.000000     0.44140        0.17466
    4     4    2.503675   -4.112800    0.000000     0.00000        0.00000
    4     5    0.489590   -0.576211    0.000000     0.11690        0.00164
    5     5   -0.695139    0.745470    0.000000     0.00000        0.00000
每一对儿激发态间的跃迁偶极矩的三个分量,彼此间的能量差,以及相应的振子强度都输出出来了。

值得一提的是,对于编号相同的输出,得到的跃迁偶极矩<m|r|n>显然就等价于这个激发态自身的偶极矩<n|r|n>。比如3 3,如上可见Multiwfn的结果为-0.434476    2.840406    0.000000,换算成Debye为单位就是-1.104329 7.219591 0.000000。Gaussian直接输出的3号激发态的偶极矩为
X=  -1.1033    Y=  7.2217    Z=  0.0000
可见相符得很好。同时这也表明IOp(9/40=3)的精度就已经挺好了。

PS:在Gaussian中输出激发态偶极矩是使用例如下面这样的关键词,在输出文件末尾就会由L601模块输出第3个激发态的偶极矩
#P b3lyp/6-31+G* TD(root=3,nstates=5) density=rhoci
注意这里用了rhoci关键词,如果只写density的话,那么Gaussian传递给L601模块的激发态的密度是弛豫的密度。而当用rhoci时,传递的是非弛豫的密度,这才是和通过组态系数和基态轨道直接产生的密度直接对应的。

评论已关闭