通过Multiwfn计算各个轨道的偶极矩

通常人们只是计算分子整体的偶极矩，有时也讨论局部的偶极矩，诸如原子偶极矩或ELF盆的偶极矩，这个靠Multiwfn (http://sobereva.com/multiwfn)的模糊空间分析和盆分析模块方便地实现。也有人问怎么计算分子轨道的偶极矩，实际上利用Multiwfn里的多个功能都可以实现这一目的，下面就依次介绍下，顺便也使读者了解Multiwfn的灵活性。例子都用B3LYP/6-31G**下的水分子。注意，由于轨道的单极矩不为0，所以轨道偶极矩的数值取决于原点的选择，下文的方法在计算的时候原点就是输入文件里的坐标原点。

(1) 最简单也是最好的方法就是利用Multiwfn 3.3.5版以上的主功能200里的选项10。这个功能能够计算轨道间的电偶极矩积分、磁偶极矩积分、速度积分、动能积分和重叠积分，详见手册3.200.10节。

启动Multiwfn后依次输入
H2O.wfn  //.wfx、.fch、.molden、NBO plot文件也都可以作为输入文件，生成它们的方法见手册第四章开头
200
10
1  //计算电偶极矩积分
4   //计算相同轨道间的积分
结果马上输出到了当前目录下的orbint.txt文件中，内容如下（目前最新版Multiwfn还有第6列，是偶极矩矢量的模）
       1       1          0.00000000      0.00000000     -0.22527032
       2       2          0.00000000     -0.00000000      0.16388523
       3       3          0.00000000      0.00000000      0.19931608
       4       4          0.00000000      0.00000000     -0.36217654
       5       5          0.00000000      0.00000000     -0.17784638
前两列是计算积分的两个轨道的编号，后三列是偶极矩XYZ分量(a.u.)。当前这个体系的轨道都是双占据的，因此将积分值乘上2就是各个轨道的偶极矩了。若把最后一列加和，乘上2，就得到了体系的总偶极矩-0.80418 a.u.。

由于.wfn和.wfx文件只记录占据轨道的信息，上例只有5个轨道。如果还要得到空轨道的偶极矩，需要用.fch或.molden文件作为输入。


(2) 利用主功能15（模糊空间分析模块）的选项2可以计算每个原子的电荷、偶极矩、四极矩、八极矩以及对分子偶极矩的贡献，并且在最后还会输出整个体系的偶极矩。而主功能6（修改波函数）的选项26可以修改轨道占据数，这会影响接下来的计算的结果。因此，若想计算某个轨道的偶极矩，就把除了这个轨道之外的其余轨道的占据数都设为0来消除它们的贡献，然后计算出的体系偶极矩扣除原子核产生的偶极矩就是这个轨道的偶极矩。

假设这里计算第4个轨道的偶极矩。启动Multiwfn后依次输入
H2O.wfn  //.wfx、.fch、.molden亦可
6  //修改波函数
26  //修改占据数
0  //选择所有轨道
0  //将占据数设为0
4  //选择第4号轨道
2  //将此轨道占据数设为它原本的数值2
q
-1  //返回主菜单
15  //模糊空间分析
2
输出最后看到体系的偶极矩
Molecular dipole moment (a.u.):       0.000000     -0.000000     -0.724353
由于Gaussian会自动把体系坐标调整到标准朝向下，原点是原子核电荷中心，所以原子核产生的偶极矩精确为0，故上面的数值就是第4号轨道的偶极矩。


(3) 主功能200里的功能2可以在Hilbert空间下计算原子和键偶极矩，这个功能也会输出体系总偶极矩。用它来计算轨道偶极矩和上例一样，也是先把轨道占据数改了就行了。注意由于此功能是基于基函数信息的，而.wfn、.wfx都不含基函数信息，因此输入文件必须用.fch或.molden。这个功能不需要像模糊空间分析模块那样做数值积分，而是解析地计算，因此速度明显更快，特别是对大体系而言。

假设按照前例只保留第4个轨道的占据数，进入这个功能后在开头会看到
Molecular nuclear dipole moment (a.u.):
 X=   -0.000000  Y=    0.000000  Z=   -0.000000  Norm=    0.000000
Molecular electron dipole moment (a.u.):
 X=    0.000000  Y=    0.000000  Z=   -0.724353  Norm=    0.724353
Molecular dipole moment (a.u.):
 X=   -0.000000  Y=    0.000000  Z=   -0.724353  Norm=    0.724353
依次是原子核产生的偶极矩、电子产生的偶极矩，以及二者的加和。这里给出的电子偶极矩即是第4号轨道的偶极矩。


(4) 由手册2.7节可见，21、22、23号用户自定义函数分别是电偶极矩XYZ分量的被积函数，比如23号用户自定义函数的定义为-z*ρ(r)，因此对它进行全空间积分，就能得到体系的偶极矩Z分量。Multiwfn的主功能100的选项4可以对任意实空间函数进行数值积分，积分精度取决于径向和角度方向积分格点数（由settings.ini里的sphpot和radpot控制，越大越精确但也越耗时），以及函数特征（越平滑积分精度越高），积分算法是Becke提出的积分DFT泛函的方法，详见《密度泛函计算中的格点积分方法》（http://sobereva.com/69）。对于-z*ρ(r)这样的比较平滑的函数，默认设定下就能积分到很高精度。

这里计算第4条轨道的偶极矩Z分量。首先将settings.ini里的iuserfunc设为23来启用第23号自定义函数，然后按照前文的(2)的例子载入输入文件并修改轨道占据数，然后退回主菜单，输入
100
4  //全空间积分实空间函数
100 //用户自定义函数
结果为-0.724353 a.u.。
 
 
以上计算轨道偶极矩的方法中(1)最简单而且快捷，(4)最繁琐。

以上方法也都可以计算其它类型的偶极矩，比如自然轨道、NTO和NLMO的偶极矩，这就看你的输入文件里存的是什么轨道了。另外，通过修改占据数，结合上面提到的功能计算原子偶极矩或多极矩，也可获知轨道对原子偶极矩或多极矩的贡献。