计算化学公社

标题: 使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩 [打印本页]
作者
Author: sobereva    时间: 2020-6-23 00:30
标题: 使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩
使用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的偶极矩矢量就是计算
(, 下载次数 Times of downloads: 176)
其中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模型描述的水环境下产生。此体系分子结构如下：
(, 下载次数 Times of downloads: 237)

启动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界面调节。

(, 下载次数 Times of downloads: 203)

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

如果之后想删掉箭头，在文本窗口运行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来让箭头长度变成偶极矩大小的两倍。

得到的图像如下

(, 下载次数 Times of downloads: 204)

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

作者
Author: yejlyejl    时间: 2020-6-23 01:32
请教一下：Multiwfn可以计算分子片段的电四极矩吗？
作者
Author: zsu007    时间: 2020-6-23 08:54
谢谢社长的分享！
作者
Author: sobereva    时间: 2020-6-23 12:29
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
作者
Author: lijiayisjtu    时间: 2020-6-23 14:18

谢谢社长的分享！又涨知识了
作者
Author: yejlyejl    时间: 2020-6-24 09:30
sobereva 发表于 2020-6-23 12:29
可以。本文已经体现了：

Molecular quadrupole moments (Standard Cartesian form):

谢谢社长指教!
再请教一下：在核磁某些原子的电四级矩效应会导致相邻核的横向驰豫时间变短，使谱峰变宽；
这个四级矩效应的大小就是这里的Molecular quadrupole moments Magnitude吗？
作者
Author: sobereva    时间: 2020-6-24 16:24
yejlyejl 发表于 2020-6-24 09:30
谢谢社长指教!
再请教一下：在核磁某些原子的电四级矩效应会导致相邻核的横向驰豫时间变短，使谱峰变宽 ...

不是
分子多极矩和核的多极矩完全是两码事
作者
Author: yejlyejl    时间: 2020-6-25 17:34
sobereva 发表于 2020-6-24 16:24
不是
分子多极矩和核的多极矩完全是两码事

谢谢！
核四极矩可以用什么软件计算吗？
作者
Author: sobereva    时间: 2020-6-27 01:13
yejlyejl 发表于 2020-6-25 17:34
谢谢！
核四极矩可以用什么软件计算吗？

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

多谢社长！
作者
Author: 哈库那摩塔塔    时间: 2022-9-1 17:17
请问软件可以分析VASP计算的二维材料中负载的团簇的偶极矩吗
作者
Author: sobereva    时间: 2022-9-2 07:38
哈库那摩塔塔 发表于 2022-9-1 17:17
请问软件可以分析VASP计算的二维材料中负载的团簇的偶极矩吗

不支持VASP
Multiwfn对周期性波函数仅支持CP2K，在可预见的未来不会支持VASP
作者
Author: 乐平    时间: 2022-11-16 08:49
社长您好！
我按照教程完成了练习，操作没问题。

但是，自己的体系里原子顺序是间断的，在用 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 作为片段，应该如何输入呢？

盼指点，谢谢！
作者
Author: sobereva    时间: 2022-11-17 01:44
乐平 发表于 2022-11-16 08:49
社长您好！
我按照教程完成了练习，操作没问题。

把-替换成[空格]to[空格]，把逗号替换成[空格]
作者
Author: 乐平    时间: 2022-11-17 11:43
本帖最后由 乐平 于 2022-11-17 11:50 编辑
sobereva 发表于 2022-11-17 01:44
把-替换成[空格]to[空格]，把逗号替换成[空格]

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



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

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



作者
Author: wanan    时间: 2023-10-19 20:06
sob老师，请问总的偶极矩可以结构偶极矩和电子偶极矩吗？如果可以的话，怎么判断总偶极矩，结构偶极矩和电子偶极矩呢？
作者
Author: sobereva    时间: 2023-10-20 08:13
wanan 发表于 2023-10-19 20:06
sob老师，请问总的偶极矩可以结构偶极矩和电子偶极矩吗？如果可以的话，怎么判断总偶极矩，结构偶极矩和电 ...

“请问总的偶极矩可以结构偶极矩和电子偶极矩吗？”   病句
作者
Author: wanan    时间: 2023-10-20 10:07
sobereva 发表于 2023-10-20 08:13
“请问总的偶极矩可以结构偶极矩和电子偶极矩吗？”   病句

不好意思老师，看了您关于18碳环在电场下的超强调控作用，您文章里面提到18碳环可以分成内平面π电子偶极矩和外平面π电子偶极矩等，所以想问下您，偶极矩可以分成结构偶极矩和电子偶极矩吗？如果可以，怎么判断呢？
作者
Author: sobereva    时间: 2023-10-21 02:49
wanan 发表于 2023-10-20 10:07
不好意思老师，看了您关于18碳环在电场下的超强调控作用，您文章里面提到18碳环可以分成内平面π电子偶极 ...

没有“结构偶极矩”这种说法
搞清楚偶极矩的定义和计算方式。取决于核坐标和电子密度分布
作者
Author: wanan    时间: 2023-10-23 09:45
sobereva 发表于 2023-10-21 02:49
没有“结构偶极矩”这种说法
搞清楚偶极矩的定义和计算方式。取决于核坐标和电子密度分布

好滴，谢谢sob老师
作者
Author: Jack    时间: 2023-11-21 23:27
请问sob老师，第一节原理部分片段偶极矩矢量计算公式有更多的说明文档或者参考文献吗？谢谢！
作者
Author: sobereva    时间: 2023-11-22 00:12
Jack 发表于 2023-11-21 23:27
请问sob老师，第一节原理部分片段偶极矩矢量计算公式有更多的说明文档或者参考文献吗？谢谢！

没什么可参考的，显而易见
如果不懂原子权重函数是什么，看Multiwfn手册3.18节的相关介绍
作者
Author: Jack    时间: 2023-11-22 00:22
sobereva 发表于 2023-11-22 00:12
没什么可参考的，显而易见
如果不懂原子权重函数是什么，看Multiwfn手册3.18节的相关介绍

好的，谢谢sob老师
作者
Author: Uus/pMeC6H4-/キ    时间: 2025-5-6 14:13
这有好几个综合性的问题：
1. 文中的计算公式涉及电子密度ρ(r)，那么在用主功能15的子功能2计算偶极矩时：
(1)如果载入multiwfn的文件已经包含格点数据，如Gaussian cube格式，此时是否直接取格点数据作为ρ(r)呢？
(2)如果载入multiwfn的文件只有波函数信息但不包含格点数据，那计算ρ(r)的中间步骤里格点数量/质量的默认设置如何呢？计算ρ(r)的空间范围是否由子选项-5定义的原子范围限制？
(3)同上，能否在计算完偶极矩并返回主菜单时把ρ(r)的格点数据保留于内存，免得后续再用主功能5算一遍？
2. 一般电子密度的格点数据是用主功能5的子功能1计算的，算完同时输出Electric dipole moment estimated by integrating electron density内容，这和主功能15的子功能2所得结果具体有什么区别？
3. 本文计算偶极矩的方法明确支持CP2K做周期性计算输出的molden文件（按http://sobereva.com/651的说法补充了盒子信息和有效核电荷数）吗？尝试载入molden后进入主功能15，发现2 Calculate atomic and molecular multipole moments and <r^2>没有显示出来（子功能1下面就是子功能3），但输入2似乎还是能用的。

问这些是因为想起来http://bbs.keinsci.com/thread-50892-1-1.html和http://bbs.keinsci.com/thread-52757-1-1.html提及的通过动力学获得红外光谱。对这种计算偶极自相关函数的目的，让程序自动输出每步偶极矩，和让程序逐帧输出波函数信息再用multiwfn批量计算，会有精度上明显的区别吗？
作者
Author: sobereva    时间: 2025-5-6 23:49
Uus/pMeC6H4-/キ 发表于 2025-5-6 14:13
这有好几个综合性的问题：
1. 文中的计算公式涉及电子密度ρ(r)，那么在用主功能15的子功能2计算偶极矩时 ...

1 (1) 不。此时选这个没意义
(2) 此功能目前用的是原子中心格点，每个原子对应的格点数是settings.ini里的radpot和sphpot设的（默认会恰当剪裁节约时间）
(3) 不能。本身用的格点形式也不一样。主功能5是均匀格点

2 主功能15的子功能2目前不支持周期性体系，这是为什么界面没显示。对于孤立体系，主功能15的子功能2是基于原子中心格点，而主功能5基于均匀格点，后者对于周期靠后的元素即便用很精细的格点也不可能算准偶极矩，因为均匀格点没法描述好它们变化剧烈的内核电子

建议直接让动力学程序输出每帧的偶极矩。




欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3