使用Multiwfn+VMD繪制片段貢獻的躍遷偶極矩矢量

1 前言

躍遷偶極矩是討論電子激發最關鍵的量之一。躍遷偶極矩大，振子強度才可能大，對應的吸收/發射才可能比較強，這些基本知識在此帖都有介紹：《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）。量化程序給的是整個分子的躍遷偶極矩，但是我們往往想考察體系中不同部分對躍遷偶極矩的貢獻，從而更好地了解躍遷的內在特征。在Multiwfn程序的主功能18里的“空穴-電子”分析模塊中，可以通過Mulliken劃分，把躍遷偶極矩分解為基函數的貢獻，以及原子的貢獻。將原子貢獻加和就可以得到某個片段的貢獻。如果我們把數據讀入到VMD程序中，還可以在VMD中繪制箭頭，直觀地考察各個片段的躍遷偶極矩矢量，使各個片段對躍遷偶極矩的貢獻一目了然。本文就介紹如何實現。
 
對Multiwfn不了解者請參考《Multiwfn FAQ》（http://www.shanxitv.org/452）、《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）。本文對應的是Multiwfn最新版本的情況，程序可以在主頁http://www.shanxitv.org/multiwfn免費下載。VMD使用的是1.9.3版，可以在http://www.ks.uiuc.edu/Research/vmd/免費下載。
 
 

2 原理

按照Mulliken劃分，第r基函數對躍遷偶極矩的貢獻寫為（注意是個矢量）
其中P_tran是某兩個態之間的躍遷密度矩陣，后面那項是r與s基函數間的電偶極矩積分。把基函數貢獻按照原子加和，就得到了原子對躍遷偶極矩的貢獻。

在《分子軌道成分的計算》（http://sioc-journal.cn/Jwk_hxxb/CN/abstract/abstract340458.shtml）一文中，筆者就已經明確強調使用Mulliken劃分時基組不能帶彌散函數，否則結果缺乏物理意義。因此，按照本文方法以Mulliken劃分來分解躍遷偶極矩時，做電子激發計算時候用的基組也絕對不能帶彌散函數。實際上，除非里德堡激發，否則帶上彌散函數對電子激發計算結果也沒明顯好處，基本是浪費時間，這里已經提過了：《亂談激發態的計算方法》（http://www.shanxitv.org/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目錄下：loaddip.tcl。然后在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正方向）

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，就得到了下圖，綠色箭頭描述體系總躍遷偶極矩，其具體數值也顯示在了文本窗口中。

要想刪除所有已繪制的箭頭，輸入draw delete all。下面我們輸入dipatm，看看各個原子貢獻的躍遷偶極矩的箭頭。由于原子貢獻的往往比較小，為了避免箭頭被原子球遮擋，drawing method我們改成line。
由圖可見，雖然每個苯環都在Y的負方向對躍遷偶極矩有巨大貢獻，但是分解到原子層面，苯環內的部分原子對躍遷偶極矩卻是沖著Y的正方向。

前面提到，片段對躍遷偶極矩的貢獻一般是有原點依賴性的。為了展示這一點，對前面的.gjf文件里每個原子Y坐標都加上6埃，并且用了nosymm關鍵詞避免自動被Gaussian平移。重復之前的作圖步驟，這次得到的圖像如下

和之前的圖對比，我們看到綠色箭頭，即總躍遷偶極矩并沒有因為體系的平移而發生改變，但是兩個苯環產生的貢獻卻與之前截然不同，這是因為這兩個苯環的躍遷電荷都不為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后面的值來調整。