使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩

sobereva

使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩

Calculation of dipole moment of molecular fragments and dipole moment of monomers in complexes using Multiwfn

文/Sobereva@北京科音  2020-Jun-23

1 原理

时不时有人问怎么计算一个分子中某些片段的偶极矩，或者问怎么计算一个分子复合物中各个分子的偶极矩，确实这对于分析局部的极性很有意义，在本文就通过例子对做法进行说明。做这种分析必须使用2020-Jun-22之后更新的Multiwfn程序，可以在http://sobereva.com/multiwfn免费下载，相关基本知识看《Multiwfn FAQ》（http://sobereva.com/452）和《Multiwfn入门tips》（http://sobereva.com/167）。此功能需要使用含有波函数信息的文件作为输入文件，比如wfn、wfx、mwfn、fch、molden等，产生的方式见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。

简单来说，计算一个体系中由某些原子构成的片段F的偶极矩矢量就是计算

其中R是原子核位置，Z是原子核电荷，ρ是电子密度，w是原子的权重函数。计算片段偶极矩矢量需要用Multiwfn的主功能15（模糊空间分析），目前支持Hirshfeld、Becke、Hirshfeld-I形式的原子权重函数。本文就用比较常用，而且物理意义比较清楚的Hirshfeld权重函数来算。

需要注意的是，如果某个片段的净电荷不为零，则这个片段的偶极矩是有原点依赖性的，即平移体系会令它发生改变。多数情况下被计算的片段的净电荷不为0，却也没有办法完全避免原点依赖性，只能是在讨论的时候留意这个问题，予以恰当讨论。

笔者之前还写过一篇与本文有一定关系的文章《使用Multiwfn+VMD绘制片段贡献的跃迁偶极矩矢量》（http://sobereva.com/396），感兴趣者推荐看看，这对于考察体系的不同区域对电子激发的振子强度的影响很有益。


2 计算分子片段的偶极矩实例

本节计算下面这个单重态[Ru(bpy)3]2+阳离子体系的各个配体的偶极矩，其wfn文件可以在这里下载：http://sobereva.com/attach/558/Ru_bpy_3.zip。此文件的结构使用BP86泛函结合SDD和6-31G*进行优化，波函数用B3LYP结合SDD和6-311G*在IEFPCM模型描述的水环境下产生。此体系分子结构如下：

启动Multiwfn后载入Ru_bpy_3.wfn，然后依次输入
15  //模糊空间分析
-1  //选择原子权重函数
3   //基于内置原子密度的Hirshfeld权重函数
-5  //定义功能1、2计算时考虑的原子范围
2-21  //体系中一个配体里的原子序号
2  //计算原子和分子多极矩

此时Multiwfn依次计算各个原子的多极矩，最后给出这些原子对应的片段的电子数、偶极矩和不同形式的多极矩：

             *****  Molecular dipole and multipole moments  *****
Total number of electrons:     81.427849
Molecular dipole moment (a.u.):      -0.000000      4.348703      0.000000
Molecular dipole moment (Debye):     -0.000000     11.053300      0.000000
Magnitude of molecular dipole moment (a.u.&Debye):      4.348703     11.053300
Molecular quadrupole moments (Standard Cartesian form):
XX=  -41.258381  XY=   -0.000000  XZ=  -15.828189
YX=   -0.000000  YY=  -13.049129  YZ=   -0.000000
ZX=  -15.828189  ZY=   -0.000000  ZZ=  -33.042376
Molecular quadrupole moments (Traceless Cartesian form):
XX=  -18.212629  XY=   -0.000000  XZ=  -23.742283
YX=   -0.000000  YY=   24.101250  YZ=   -0.000000
ZX=  -23.742283  ZY=   -0.000000  ZZ=   -5.888621
Magnitude of the traceless quadrupole moment tensor:   25.129610
Molecular quadrupole moments (Spherical harmonic form):
Q_2,0 =  -5.888621   Q_2,-1=  -0.000000   Q_2,1= -27.415227
Q_2,-2=  -0.000000   Q_2,2 = -24.429930
Magnitude: |Q_2|=   37.189945
Molecular octopole moments (Cartesian form):
XXX=    0.0000  YYY= -455.8273  ZZZ=   -0.0000  XYY=   -0.0000  XXY= -217.0765
XXZ=    0.0000  XZZ=    0.0000  YZZ= -176.7460  YYZ=    0.0000  XYZ=  -85.0737
Molecular octopole moments (Spherical harmonic form):
Q_3,0 =    -0.0000  Q_3,-1=   -20.8698  Q_3,1 =     0.0000
Q_3,-2=  -329.4892  Q_3,2 =    -0.0000  Q_3,-3=  -154.4790  Q_3,3 =     0.0000
Magnitude: |Q_3|=    364.5030

这些量的具体定义看Multiwfn手册3.18.3节，这里不细说。我们当前关心的只有
Molecular dipole moment (a.u.):      -0.000000      4.348703      0.000000
这一行，即偶极矩矢量为(0.0,4.3487,0.0) a.u.。

然后我们计算其它片段的偶极矩，依次输入
-5  //重新定义片段
42-61
2  //计算原子和分子多极矩
-5  //重新定义片段
22-41
2  //计算原子和分子多极矩

汇总一下，三个配体的偶极矩分别为：
2-21片段： ( 0.000000  4.348703  0.000000) a.u.
42-61片段：( 3.766463 -2.173911  0.001788) a.u.
22-41片段：(-3.766463 -2.173911 -0.001788) a.u.

下面，我们用VMD程序将偶极矩的箭头画出来以便于观看。VMD可以在http://www.ks.uiuc.edu/Research/vmd/免费下载，本文用的是VMD 1.9.3。wfn格式没法载入VMD，所以我们先把当前体系转换为VMD可以载入的xyz格式。在当前Multiwfn窗口中输入
0  //返回主菜单
100  //其它功能(Part 1)
2   //导出文件
2   //导出为xyz格式
直接按回车用默认路径，然后Ru_bpy_3.xyz就产生在当前目录下了。

将Ru_bpy_3.xyz载入VMD，然后将Multiwfn的examples\scripts目录下的drawarrow.tcl里的全部内容复制到VMD文本窗口里，这样就定义了一个新的名为drawarrow的命令，我们用这个命令来绘制片段偶极矩。

在VMD的文本窗口里运行以下命令
draw color green
drawarrow "serial 2 to 21" 0.000000 4.348703 0.000000 0.7
draw color red
drawarrow "serial 42 to 61" 3.766463 -2.173911  0.001788 0.7
draw color yellow
drawarrow "serial 22 to 41" -3.766463 -2.173911 -0.001788 0.7

这些命令代表分别用绿、红、黄色箭头把三个片段偶极矩绘制出来。双引号里面是选择相应片段里的原子的选择语句，语法详见《VMD里原子选择语句的语法和例子》（http://sobereva.com/504）。之后的三个参数是偶极矩的X/Y/Z分量，最后的参数0.7是在绘图时把偶极矩数值乘上这个系数再绘制，起到控制箭头长度、令图像美观的目的。绘制出来的图像如下。分子的显示方式可以在Graphics - Representation界面调节。

偶极矩是从负电中心指向正电中心。由图可见每个片段的始端都是两个氮的方向，这是因为氮带负电，而且有丰富的孤对电子。显然当前的结果完全合理。

如果之后想删掉箭头，在文本窗口运行draw delete all命令。

值得一提的是，当前这个体系里三个片段的偶极矩是对称分布的，这是因为中心金属正好在坐标原点上。由于每个配体的净电荷不为零，因此如果在做产生波函数的任务之前把体系平移一下，最终得到的三个箭头就不以Ru为中心对称了。所以要注意原点位置问题。
有读者问如果片段里的原子序号是不连续的怎么输入。假设片段里原子序号是4,9-15,89-100,111，上面例子里的双引号里就写serial 4 9 to 15 89 to 100 111，即把-替换为[空格]to[空格]，把逗号替换为[空格]，即可。


3 计算分子复合物中单体偶极矩实例

下面这个例子我们计算一下苯酚二聚体中每个苯酚单体的偶极矩，并且连同二聚体的总偶极矩一起用VMD画出来。

启动Multiwfn，然后输入
examples\phenoldimer.wfn  //自带的苯酚二聚体波函数文件
15  //模糊空间分析
-1  //选择原子权重函数
3   //Hirshfeld
2  //计算原子和分子多极矩
由于我们当前还没定义片段，所以Multiwfn最后给出的下面的信息直接就是二聚体的偶极矩了
Molecular dipole moment (a.u.):       1.227306     -0.128087      0.650833

之后和上一节的例子一样，对第一个苯酚（1~13号原子）和第二个苯酚（14~26号原子）分别计算片段偶极矩，结果分别为(0.570356 -0.356257 0.284025) a.u.和(0.656950 0.228171 0.366808) a.u.。这两个苯酚的偶极矩的矢量和与前面输出的二聚体的偶极矩是相同的。

和前例一样，把当前结构导出成xyz文件并载入VMD，并且把drawarrow.tcl里的内容拷到VMD的文本窗口里运行，之后输入以下命令进行绘制：
draw color green
drawarrow all 1.227306 -0.128087 0.650833 2
draw color red
drawarrow "serial 1 to 13" 0.570356 -0.356257 0.284025 2
draw color yellow
drawarrow "serial 14 to 26" 0.656950 0.228171 0.366808 2

这些语句代表二聚体的偶极矩用绿色箭头绘制，all代表整个体系。两个苯酚的偶极矩分别用红色和黄色绘制。由于苯酚的偶极矩较小，为了让图中箭头足够长，所以命令末尾用了2来让箭头长度变成偶极矩大小的两倍。

得到的图像如下

可见由于两个苯酚的偶极矩的方向基本相同，二者的矢量和所对应的二聚体的偶极矩比单体偶极矩大得多。

yejlyejl

请教一下：Multiwfn可以计算分子片段的电四极矩吗？

zsu007

谢谢社长的分享！

sobereva

yejlyejl 发表于 2020-6-23 01:32
请教一下：Multiwfn可以计算分子片段的电四极矩吗？

可以。本文已经体现了：

Molecular quadrupole moments (Standard Cartesian form):
XX=  -41.258381  XY=   -0.000000  XZ=  -15.828189
YX=   -0.000000  YY=  -13.049129  YZ=   -0.000000
ZX=  -15.828189  ZY=   -0.000000  ZZ=  -33.042376
Molecular quadrupole moments (Traceless Cartesian form):
XX=  -18.212629  XY=   -0.000000  XZ=  -23.742283
YX=   -0.000000  YY=   24.101250  YZ=   -0.000000
ZX=  -23.742283  ZY=   -0.000000  ZZ=   -5.888621
Magnitude of the traceless quadrupole moment tensor:   25.129610
Molecular quadrupole moments (Spherical harmonic form):
Q_2,0 =  -5.888621   Q_2,-1=  -0.000000   Q_2,1= -27.415227
Q_2,-2=  -0.000000   Q_2,2 = -24.429930
Magnitude: |Q_2|=   37.189945

lijiayisjtu

谢谢社长的分享！又涨知识了

yejlyejl

sobereva 发表于 2020-6-23 12:29
可以。本文已经体现了：

Molecular quadrupole moments (Standard Cartesian form):

谢谢社长指教!
再请教一下：在核磁某些原子的电四级矩效应会导致相邻核的横向驰豫时间变短，使谱峰变宽；
这个四级矩效应的大小就是这里的Molecular quadrupole moments Magnitude吗？

sobereva

yejlyejl 发表于 2020-6-24 09:30
谢谢社长指教!
再请教一下：在核磁某些原子的电四级矩效应会导致相邻核的横向驰豫时间变短，使谱峰变宽 ...

不是
分子多极矩和核的多极矩完全是两码事

yejlyejl

sobereva 发表于 2020-6-24 16:24
不是
分子多极矩和核的多极矩完全是两码事

谢谢！
核四极矩可以用什么软件计算吗？

sobereva

yejlyejl 发表于 2020-6-25 17:34
谢谢！
核四极矩可以用什么软件计算吗？

我不是很清楚，这应该是可以直接查到的。如果真是需要算的话，那是核物理方面的程序了，和量化没有关系

yejlyejl

sobereva 发表于 2020-6-27 01:13
我不是很清楚，这应该是可以直接查到的。如果真是需要算的话，那是核物理方面的程序了，和量化没有关系

多谢社长！

哈库那摩塔塔

请问软件可以分析VASP计算的二维材料中负载的团簇的偶极矩吗

sobereva

哈库那摩塔塔 发表于 2022-9-1 17:17
请问软件可以分析VASP计算的二维材料中负载的团簇的偶极矩吗

不支持VASP
Multiwfn对周期性波函数仅支持CP2K，在可预见的未来不会支持VASP

乐平

社长您好！
我按照教程完成了练习，操作没问题。

但是，自己的体系里原子顺序是间断的，在用 VMD 绘制箭头的时候出现报错。

1. vmd > source drawarrow.tcl
   
2. vmd > draw color green
   
3. 0
   
4. vmd > drawarrow all -0.000002     -0.000019      0.000028
   
5. Geometry center: -6.905773625476286e-5 0.001978856511414051 -0.0008773042936809361
   
6. 2
   
7. vmd > draw color red
   
8. 3
   
9. vmd > drawarrow "serial 1-36,73-74,87-94,103-112" -0.021702     -0.000006      0.000864
   
10. ERROR) Selection terminated too early
    
11. ERROR) syntax error
    
12. ERROR) Bad character:44:,
    
13. ERROR) Bad character:44:,
    
14. ERROR) Bad character:44:,
    
15. atomselect: cannot parse selection text: serial 1-36,73-74,87-94,103-112

复制代码

也尝试了如下输入方式：

1. vmd > drawarrow "serial 1 to 36, serial 73 to 74, serial 87 to 94, serial 103 to 112" -0.021702     -0.000006      0.000864
   
2. ERROR) Bad character:44:,
   
3. ERROR) Selection terminated too early
   
4. ERROR) syntax error
   
5. ERROR) Bad character:44:,
   
6. ERROR) Bad character:44:,
   
7. atomselect: cannot parse selection text: serial 1 to 36, serial 73 to 74, serial 87 to 94, serial 103 to 112
   
8. vmd >
   
9. vmd >

复制代码

依旧报错……


是不是 drawarrow.tcl 代码只支持一整段连续编号的片段？
比如，您的例子中的

1. drawarrow "serial 1 to 13" 0.570356 -0.356257 0.284025 2

复制代码

如果要将我的体系中 1-36,73-74,87-94,103-112 作为片段，应该如何输入呢？

盼指点，谢谢！

sobereva

乐平 发表于 2022-11-16 08:49
社长您好！
我按照教程完成了练习，操作没问题。

把-替换成[空格]to[空格]，把逗号替换成[空格]

乐平

sobereva 发表于 2022-11-17 01:44
把-替换成[空格]to[空格]，把逗号替换成[空格]

感谢社长指点，可以绘制了。



顺带发现了删除绘制错误箭头的方法：
draw delete all  # 删除绘制的所有箭头

draw delete 数字编号 # 删除数字编号对应的箭头（因为在输入 draw 绘制命令后回车，会显示对应的数字编号，根据数字编号可以删除对应的箭头）