24小时热门版块排行榜

返回列表

【奖励】本帖被评价4次，作者sobereva增加金币 3.2 个

sobereva

至尊木虫 (著名写手)

[资源] 利用Multiwfn计算Gaussian输出的激发态之间的跃迁偶极矩已有1人参与

利用Multiwfn计算Gaussian输出的激发态之间的跃迁偶极矩
文/Sobereva(3)
First release: 2014-Mar-28 Last update: 2014-Apr-14

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

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

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就够了，精度足够，计算速度也比较快。x=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)的精度足够了，没太大必要用到IOp(9/40=4)。

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

[ Last edited by sobereva on 2014-4-14 at 19:08 ]

回复此楼