计算化学公社
标题: 使用Multiwfn+VMD绘制片段贡献的跃迁偶极矩矢量 [打印本页]
作者Author: sobereva 时间: 2017-12-3 00:14
标题: 使用Multiwfn+VMD绘制片段贡献的跃迁偶极矩矢量
使用Multiwfn+VMD绘制片段贡献的跃迁偶极矩矢量
Using Multiwfn+VMD to plot transition dipole moment vector contributed by specific fragment
文/Sobereva @北京科音 First release: 2017-Dec-2 Last update: 2020-Dec-25
1 前言
跃迁偶极矩是讨论电子激发最关键的量之一。跃迁偶极矩大，振子强度才可能大，对应的吸收/发射才可能比较强，这些基本知识在此帖都有介绍：《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（http://sobereva.com/314）。量化程序给的是整个分子的跃迁偶极矩，但是我们往往想考察体系中不同部分对跃迁偶极矩的贡献，从而更好地了解跃迁的内在特征。在Multiwfn程序的主功能18里的“空穴-电子”分析模块中，可以通过Mulliken划分，把跃迁偶极矩分解为基函数的贡献，以及原子的贡献。将原子贡献加和就可以得到某个片段的贡献。如果我们把数据读入到VMD程序中，还可以在VMD中绘制箭头，直观地考察各个片段的跃迁偶极矩矢量，使各个片段对跃迁偶极矩的贡献一目了然。本文就介绍如何实现。
对Multiwfn不了解者请参考《Multiwfn FAQ》（http://sobereva.com/452）、《Multiwfn入门tips》（http://sobereva.com/167）和《Multiwfn波函数分析程序的意义、功能与用途》（http://sobereva.com/184）。本文对应的是Multiwfn最新版本的情况，程序可以在主页http://sobereva.com/multiwfn免费下载。VMD使用的是1.9.3版，可以在http://www.ks.uiuc.edu/Research/vmd/免费下载。
2 原理
按照Mulliken划分，第r基函数对跃迁偶极矩的贡献写为（注意是个矢量）
(, 下载次数 Times of downloads: 123)
其中P_tran是某两个态之间的跃迁密度矩阵，后面那项是r与s基函数间的电偶极矩积分。把基函数贡献按照原子加和，就得到了原子对跃迁偶极矩的贡献。
在《分子轨道成分的计算》（http://sioc-journal.cn/Jwk_hxxb/CN/abstract/abstract340458.shtml）一文中，笔者就已经明确强调使用Mulliken划分时基组不能带弥散函数，否则结果缺乏物理意义。因此，按照本文方法以Mulliken划分来分解跃迁偶极矩时，做电子激发计算时候用的基组也绝对不能带弥散函数。实际上，除非里德堡激发，否则带上弥散函数对电子激发计算结果也没明显好处，基本是浪费时间，这里已经提过了：《乱谈激发态的计算方法》（http://sobereva.com/265）。
要注意，虽然体系总的跃迁偶极矩没有原点依赖性，但是片段对跃迁偶极矩的贡献往往是有原点依赖性的。换句话说，在计算时把体系进行平移，不影响总跃迁偶极矩的结果，但是可能会影响片段产生的贡献。然而原点应设在哪里，完全是任意的，因此在讨论时一定要加以注意这个问题（Gaussian计算时如果没用nosymm关键词，默认会平移体系使原点在体系原子核电荷中心位置）。之所以会有这个特征，是因为容易证明，只有最低阶的非零多极矩才没有原点依赖性，比如中性体系偶极矩没有原点依赖性，但是离子体系由于单极矩不为0，因此偶极矩有原点依赖性，因此没太大意义。类似地，对于电子跃迁问题，体系总的跃迁电荷为0，即∑i∑j P_tran(i,j)*S(i,j)=0，因此总跃迁偶极矩没有原点依赖性，而某个基函数i的跃迁电荷往往不为零，即∑j P_tran(i,j)*S(i,j)≠0，因此基函数、原子、分子片段对跃迁偶极矩的贡献往往是有原点依赖性的。
3 实例：偶氮苯(Azobenzene)
这里以一个简单分子偶氮苯为例，介绍一下将Multiwfn与VMD相结合，绘制各个片段跃迁偶极矩矢量的完整流程。
以下是偶氮苯在PBE0/6-31G*级别做TDDFT电子激发计算的Gaussian输入文件，计算最低5个激发态。其中9/40=4必须加，否则只有较大的组态系数才会输出出来，会导致之后Multiwfn产生的跃迁偶极矩不准确（要稍微更准确的结果可以写9/40=5，不过输出文件会更大）。
%chk=C:\Azobenzene.chk
# pbe1pbe/6-31g(d) TD(nstates=5) iop(9/40=4)
pbe1pbe/6-31g(d) opted
0 1
C -0.18984100 4.51585300 0.00000000
C 1.10113800 3.99491100 0.00000000
C 1.29053000 2.61740200 0.00000000
C 0.18984100 1.75735600 0.00000000
C -1.10944900 2.28225900 0.00000000
C -1.29157500 3.65642100 0.00000000
H -0.34218800 5.59191200 0.00000000
H 1.95936300 4.66120300 0.00000000
H 2.28423900 2.17907300 0.00000000
H -1.94924600 1.59527000 0.00000000
H -2.29790100 4.06725100 0.00000000
N 0.49896400 0.37896400 0.00000000
N -0.49896400 -0.37896400 0.00000000
C -0.18984100 -1.75735600 0.00000000
C -1.29053000 -2.61740200 0.00000000
C 1.10944900 -2.28225900 0.00000000
C -1.10113800 -3.99491100 0.00000000
H -2.28423900 -2.17907300 0.00000000
C 1.29157500 -3.65642100 0.00000000
H 1.94924600 -1.59527000 0.00000000
C 0.18984100 -4.51585300 0.00000000
H -1.95936300 -4.66120300 0.00000000
H 2.29790100 -4.06725100 0.00000000
H 0.34218800 -5.59191200 0.00000000
将Azobenzene.chk用formchk转换为Azobenzene.fch。然后启动Multiwfn，载入Azobenzene.fch，之后依次输入：
18 //电子激发分析功能
11 //将跃迁偶极矩分解为基函数和原子的贡献
Azobenzene.out //Gaussian输出文件
2 //要考察的激发态。我们随便选一个，比如第2激发态
1 //分解的是跃迁电偶极矩
n //不生成AAtrdip.txt（原子-原子跃迁偶极矩矩阵）
当前目录下得到了trdipcontri.txt，其中第一部分是各个基函数对跃迁偶极矩X,Y,Z分量的贡献，第二部分是各个原子的贡献。将trdipcontri.txt拷到VMD目录下。
退到程序主菜单，然后进入主功能100的子功能2，选择导出体系的pdb文件。启动VMD，将体系的pdb文件拖入VMD main窗口载入。
将这个loadip.tcl拷到VMD目录下：
(, 下载次数 Times of downloads: 205)
。然后在VMD的文本窗口运行source loaddip.tcl命令执行之。这会将VMD目录下的trdipcontri.txt里记录的各个原子对跃迁偶极矩的贡献载入到内存（同时也输出到了文本窗口里），还定义了dip命令，用来绘制箭头表现某个片段的跃迁偶极矩，还同时定义了dipatm命令，可以一次性绘制每个原子的跃迁偶极矩箭头。
下面，我们把偶氮苯第一个苯环（原子序号1~11）、中间两个氮（原子序号12、13）和另一个苯环（原子序号14~24）分别作为三个片段来绘制它们贡献的跃迁偶极矩。用dip命令默认绘制的箭头是蓝色的，为了更显眼，先在VMD窗口里输入draw color red来让之后绘制的物体成为红色，然后在VMD文本窗口中依次输入
dip "serial 1 to 11"
dip "serial 12 13"
dip "serial 14 to 24"
这里诸如serial 1 to 11代表选择1到11号原子。dip命令绘制的箭头的圆柱部分长度对应片段贡献的跃迁偶极矩大小，箭头方向表示跃迁偶极矩矢量方向，箭头的中央位置对应选中的片段的几何中心。每次用dip命令后，在VMD的文本窗口中还会把这个片段的几何中心以及对跃迁偶极矩贡献的X,Y,Z分量直接输出出来。
注：被绘制的片段里的原子序号可以不连续。例如运行dip "serial 1 5 to 8 11 to 14 18"，会对由1、5、6、7、8、11、12、13、14、18号原子构成的片段进行绘制。
在文本窗口输入color Display Background white把背景设成白色，然后进graphics-representation，把Drawing method设为CPK。笔者建议用正交视角观看当前图像，以免因近大远小造成视觉误差，故选Display-Orthographic。最终看到的效果如下所示（左下角的坐标轴的红、绿、蓝分别对应X,Y,Z正方向）
(, 下载次数 Times of downloads: 133)
Gaussian输出的此体系S0->S2的总跃迁偶极矩，以及VMD文本窗口中看到的苯环1、N2、苯环2的贡献分别为（a.u.）
Total：0.1155 -2.8868 0.0
苯环1：0.14262 -1.43288 0.0
N2： -0.16948 -0.02138 0.0
苯环2：0.14262 -1.43288 0.0
可见跃迁偶极矩主要是Y轴负方向的。将定量数据和图像结合观察，可见两个苯环起到了最主要的贡献。而N2对Y方向贡献甚微，在X方向有一点微小贡献。
我们也可以同时把体系总的跃迁偶极矩箭头绘制出来。为了令其颜色与片段的区分，我们先输入draw color green，然后再输入dip all，就得到了下图，绿色箭头描述体系总跃迁偶极矩，其具体数值也显示在了文本窗口中。
(, 下载次数 Times of downloads: 123)
要想删除所有已绘制的箭头，输入draw delete all。下面我们输入dipatm，看看各个原子贡献的跃迁偶极矩的箭头。由于原子贡献的往往比较小，为了避免箭头被原子球遮挡，drawing method我们改成line。
(, 下载次数 Times of downloads: 135)
由图可见，虽然每个苯环都在Y的负方向对跃迁偶极矩有巨大贡献，但是分解到原子层面，苯环内的部分原子对跃迁偶极矩却是冲着Y的正方向。
前面提到，片段对跃迁偶极矩的贡献一般是有原点依赖性的。为了展示这一点，对前面的.gjf文件里每个原子Y坐标都加上6埃，并且用了nosymm关键词避免自动被Gaussian平移。重复之前的作图步骤，这次得到的图像如下
(, 下载次数 Times of downloads: 158)
和之前的图对比，我们看到绿色箭头，即总跃迁偶极矩并没有因为体系的平移而发生改变，但是两个苯环产生的贡献却与之前截然不同，这是因为这两个苯环的跃迁电荷都不为0（用Multiwfn的空穴-电子分析界面的子功能6可以输出每个原子的跃迁电荷，按照片段进行加和，得到的两个苯环的数值分别为0.2116和-0.2116）。而N2部分跃迁电荷恰为0，因此它贡献的跃迁偶极矩和平移前一样，依然是-0.16948 -0.02138 0.0。
4 其它
上面的例子是将基态与激发态之间的跃迁电偶极矩分解成原子/片段的贡献，也可以将跃迁磁偶极矩以相同的方法进行分解并绘制，只要在Multiwfn问你分解哪种跃迁偶极矩的时候选Magnetic即可。另外，2020-Dec-25及之后更新的Multiwfn也支持对两个激发态之间分解它们的跃迁电偶极矩，让你输入激发态序号的那一步输入两个相应的激发态序号即可。
顺带一提，如果想令VMD中绘制的箭头长一点或者短一点，可以修改loaddip.tcl当中的这一部分：
set begx [expr $cenx-$fragdx/2]
set begy [expr $ceny-$fragdy/2]
set begz [expr $cenz-$fragdz/2]
set endx [expr $cenx+$fragdx/2]
set endy [expr $ceny+$fragdy/2]
set endz [expr $cenz+$fragdz/2]
将$fragdz改为比如2*$fragdz，就可以让箭头长度变为之前两倍。箭头的粗细可以通过修改draw cylinder这一句的radius后面的值来调整。
作者Author: Jasminer 时间: 2017-12-3 20:44
弱弱地问：既然跃迁偶极矩的片段划分不具有守恒性，做了有什么意义呢？
作者Author: muxijiao 时间: 2017-12-3 22:09
社长拿Azobenzene做例子，好亲切
作者Author: sobereva 时间: 2017-12-4 00:58
虽然有原点依赖性，但是按照习俗，原点放到体系中心时来讨论，还是能讨论很多问题的。
原点依赖性来自于偶极矩积分部分，有任意性，但跃迁密度矩阵是没有原点依赖性的，是严格的。
作者Author: lindlar 时间: 2017-12-5 10:48
学习啦
作者Author: Realzlzl 时间: 2017-12-12 17:26
sob老师，可能是我操作有误，按照您的操作流程，
18 //电子激发分析功能
1 //空穴-电子分析
Azobenzene.out //Gaussian输出文件
2 //要考察的激发态。我们随便选一个，比如第二激发态
5 //将跃迁偶极矩进行分解
1 //程序可以分解跃迁电偶极矩和跃迁磁偶极矩，这里选择分解前者
当前目录下得到了trdipcontri.txt，其中第一部分是各个基函数对跃迁偶极矩X,Y,Z分量的贡献，第二部分是各个原子的贡献。将trdipcontri.txt拷到VMD目录下。
选0退到程序主菜单（这儿选择0后，会出现一系列数值），然后进入主功能100的子功能2（不能进入主功能100），选择导出体系的pdb文件。启动VMD，将体系的pdb文件拖入VMD main窗口载入。
麻烦sob老师了！！！
作者Author: sobereva 时间: 2017-12-12 19:25
再输入一次0，从主功能18的菜单退回到主菜单，再进入主功能100
作者Author: Realzlzl 时间: 2017-12-13 09:22
sob老师，您好，某种分子的极化率是常数嘛？与电场或者波长有一定的函数关系吗？谢谢！
作者Author: sobereva 时间: 2017-12-13 13:59
极化率分静态和含频极化率，后者和外场频率有关，前者是0频极限的值。
作者Author: Realzlzl 时间: 2017-12-14 16:13
好的，谢谢！！！
作者Author: sobereva 时间: 2020-12-25 06:59
今日在官网更新的Multiwfn已经支持对两个激发态之间计算原子/片段对其跃迁电偶极矩的贡献
作者Author: Sway 时间: 2021-6-22 21:25
sob:老师怎样把电/磁偶极矩矢量在同一个图中表现出来?
作者Author: sobereva 时间: 2021-6-22 23:48
分清楚磁偶极矩和跃迁磁偶极矩
要在图中表现什么矢量，就用VMD中的绘图命令相应地画箭头就完了
作者Author: 桌坛新星10500 时间: 2022-2-10 21:11
"然后在VMD文本窗口中依次输入
dip "serial 1 to 11"
dip "serial 12 13"
dip "serial 14 to 24"
这里诸如serial 1 to 11代表选择1到11号原子."
sob老师您好，请教老师在这一步中，如果我要考虑的某个片段原子序号不是连续的例如：“2,5,7,12-16,19-22,27-59”为一个取代基片段的原子序号，这种情况该如何在vmd中输入呢？
作者Author: sobereva 时间: 2022-2-11 01:23
-用 to 替换，逗号用空格替换
作者Author: 桌坛新星10500 时间: 2022-2-11 23:14
谢谢社长解答！
作者Author: 静静子 时间: 2022-7-21 12:14
老师您好，按照您的帖子流程，
18 //电子激发分析功能
1 //空穴-电子分析
Azobenzene.out //Gaussian输出文件
1 //选第一激发态
5 //将跃迁偶极矩进行分解
1 //跃迁电偶极矩
当前目录下得到了trdipcontri.txt，将trdipcontri.txt拷到VMD目录下，并将帖子中的loaddip.tcl文件拷入VMD文件夹下，启动VMD，将pdb文件导入，然后在VMD的文本窗口运行source loaddip.tcl命令执行，结果显示couldn't execute "C:\Users\wenjing\Desktop\loaddip.tcl": no such file or directory，可能是我自己的操作问题，但不知是哪里出了问题，麻烦sob老师了
作者Author: sobereva 时间: 2022-7-21 17:59
VMD启动后文本窗口输入pwd，显示的目录在哪，tcl脚本就放到哪
作者Author: 静静子 时间: 2022-10-6 14:32
set begx [expr $cenx-$fragdx/2]
set begy [expr $ceny-$fragdy/2]
set begz [expr $cenz-$fragdz/2]
set endx [expr $cenx+$fragdx/2]
set endy [expr $ceny+$fragdy/2]
set endz [expr $cenz+$fragdz/2]
将$fragdz改为比如2*$fragdz，就可以让箭头长度变为之前两倍。是加减号后面的$fragdz都改为2*$fragdz，才表示键长变为两倍嘛？
作者Author: 静静子 时间: 2022-10-6 14:32
set begx [expr $cenx-$fragdx/2]
set begy [expr $ceny-$fragdy/2]
set begz [expr $cenz-$fragdz/2]
set endx [expr $cenx+$fragdx/2]
set endy [expr $ceny+$fragdy/2]
set endz [expr $cenz+$fragdz/2]
将$fragdz改为比如2*$fragdz，就可以让箭头长度变为之前两倍。是加减号后面的$fragdz都改为2*$fragdz，才表示键长变为两倍嘛？
作者Author: snljty2 时间: 2022-10-6 19:17
静静子 发表于 2022-10-6 14:32
set begx [expr $cenx-$fragdx/2]
set begy [expr $ceny-$fragdy/2]
set begz [expr $cenz-$fragdz/2]
把/2删了。
作者Author: 静静子 时间: 2022-10-13 16:44
好的，感谢
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |