通過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的偶極矩，這就看你的輸入文件里存的是什么軌道了。另外，通過修改占據數，結合上面提到的功能計算原子偶極矩或多極矩，也可獲知軌道對原子偶極矩或多極矩的貢獻。