使用Multiwfn分析Gaussian的極化率、超極化率的輸出
使用Multiwfn分析Gaussian的極化率、超極化率的輸出
文/Sobereva @北京科音
First release: 2014-Apr-27 Last update: 2021-May-14
1 前言
Gaussian中polar關鍵詞是專門用來計算極化率α和第一超極化率β的，從Gaussian09 D.01開始也支持polar=gamma計算第二超極化率γ。Gaussian輸出的很多信息總是令初學者很難看懂，（超）極化率的輸出更是如此，內容雜亂，而且不同理論方法下輸出的信息格式還差異很大，分析起來很不方便。Multiwfn的主功能200里的子功能7就是專門用來解析Gaussian的（超）極化率的輸出的，使之簡潔易懂。與此同時，Multiwfn還會計算一些對于分析（超）極化率十分重要的量，比如<α>、α的各向異性值、β在分子偶極矩方向的分量等等，這使得通過Gaussian研究（超）極化率變得方便得多。此文的目的就是介紹如何使用Multiwfn分析Gaussian的（超）極化率的輸出。
本文內容對應的是Multiwfn官網上的最新版本，可以在其主頁http://www.shanxitv.org/multiwfn上免費下載，不要用老版本。Multiwfn是最為強大的波函數分析程序，此文介紹的只不過是它的一個附屬功能而已。如果對Multiwfn不熟悉，建議閱讀《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。Multiwfn分析polar關鍵詞輸出的這個功能只支持Gaussian 09及之后版本，不支持更老版本。
Multiwfn的這個功能還支持計算與超瑞利散射(Hyper-Rayleigh Scattering, HRS)實驗相關的各種量，在另一篇文章《使用Multiwfn計算與超瑞利散射(HRS)實驗相關的量》（http://www.shanxitv.org/499）當中做了詳細介紹。
如果你使用了本文介紹的Multiwfn的功能，請在發表文章時提及研究中使用了Multiwfn提取Gaussian輸出文件中的相關數據并計算了相關的量，并應當根據Multiwfn程序啟動時提示的信息恰當引用Multiwfn程序。
2 Gaussian中（超）極化率的計算
這一節介紹在Gaussian中計算（超）極化率的基本原理和所需要的關鍵詞。
將能量E對均勻外電場F進行Taylor展開得到如下式子
μ0就是分子的永久偶極矩矢量（無外場下的偶極矩）。α是分子的極化率，也稱線性光學系數，是個3*3二階張量。β是第一超極化率，也稱分子的二階非線性光學(NLO)系數，是個3*3*3的三階張量。γ是第二超極化率，亦稱分子的三階NLO系數，研究得相對少一些。而更高的δ極少被討論。（超）極化率數值和電場頻率是相關的，外場頻率為0的情況被稱為靜態（超）極化率；如果外場頻率不是0，比如某一頻率的光，對應的就是動態（超）極化率，或者叫含頻（超）極化率。
Gaussian的polar關鍵詞求這些量是通過導數方法，其原理很明確。如上式可見，求極化率要令能量對外場做兩次導數，求第一超極化率要做三次導數。
求能量的導數分為三種情況：
（1）解析導數。這種方法求導數又快又精確，但是編程實現很困難，特別是對于高級的后HF而言。這種方式求能量導數時需要用到分子軌道系數對外場的導數，這通過求解CPHF方程來實現。并且利用含頻CPHF方程，可以得到動態（超）極化率。這種求（超）極化率的方式也被稱為CPHF方法。
（2）數值導數。這種方式通過有限差分方式得到能量的導數，也稱為有限場(FF)方法。有限差分在所有數值算法書里都會介紹，比如求一個函數f在0處的導數可表示為f(0)=[f(0.001)-f(-0.001)]/(2*0.001)。通過這種方法求（超）極化率非常費時間，需要在不同數值、不同方向的外場下計算很多很多次單點，而且求高階導數精度很差，因為每求一次數值導數都會由于數值噪音而累積誤差，隨著求導階數增加誤差迅速放大，做二階以上導數時精度就很爛了（除非通過額外修正、外推等特殊方式處理）。另外，通過數值導數也沒法求得動態（超）極化率。
(3) 解析+數值混合導數。這就是基于低階解析導數，再通過一次或多次有限差分求高階導數以獲得（超）極化率的方法。精度和速度介于全解析導數和全數值導數之間。
Gaussian在支持解析導數方面已經算是很好的，對于HF、DFT、半經驗方法都能支持到三階解析導數，因此可以全解析地得到靜態和含頻的α和β（而其它量化程序都很難對HF/DFT/半經驗支持到這么高階導數，通常也就能支持到二階）。對這類支持三階解析導數的方法，使用polar關鍵詞就會直接給出α和β的結果。如果同時寫上CPHF=RdFreq，并且在輸入文件末尾空一行寫上外場頻率w，比如0.05 0.12 0.3或532nm 670nm（靜態的情況總是會被計算，故不需要特意寫），就會計算含頻極化率α(-ω;ω)和含頻第一超極化率β(-ω;ω,0)。如果要算的含頻第一超極化率是SHG形式，即β(-2ω;ω,ω)，就寫polar=DCSHG即可，此時CPHF=RdFreq可以忽略不寫，并且此時β(-ω;ω,0)也會照樣計算。
對于MP2等Gaussian只能支持到二階解析導數的方法，可以全解析地計算α，但是若想計算β，就得再做一次有限差分才行。Gaussian中，此時直接寫polar只會計算α，若還想得到β就得寫polar=Cubic，會自動基于解析二階導數做有限差分來得到三階導數。
對于CISD、QCISD、CCSD、MP3、MP4(SDQ)等Gaussian只支持一階解析導數的方法，直接寫polar關鍵詞會基于解析一階導數做一次有限差分得到二階導數，即α。若還想得到β就得寫polar=DoubleNumer（等價于polar=enonly），會基于解析一階導數做兩次有限差分得到三階導數。
對于CCSD(T)、QCISD(T)、MP4(SDTQ)、MP5等Gaussian只支持計算單點能的方法，寫polar關鍵詞會對能量做兩次有限差分并由此得到α。對于它們在Gaussian里沒有辦法直接給出β。
注意對于HF、DFT、半經驗以外的方法，Gaussian都沒法給出含頻（超）極化率。
若想得到更高階導數計算更高階超極化率，比如DFT下計算第二超極化率γ，需要基于其解析三階導數做一次有限差分。從G09 D.01開始，對于支持三階解析導數的方法可以用polar=gamma來算靜態γ和動態的γ(-2ω;ω,ω,0)、γ(-ω;ω,0,0)，會基于解析三階導數自動做有限差分來得到所需的四階導數。末尾應空一行寫上外場頻率。
值得一提的是，利用Multiwfn，可以通過完全態求和方法基于Gaussian的CIS/TDHF/TDDFT計算的輸出得到靜態或動態的極化率和第一、二、三超極化率，見《使用Multiwfn基于完全態求和(SOS)方法計算極化率和超極化率》（http://www.shanxitv.org/232）。Multiwfn還能圖形化展現（超）極化率，見《使用Multiwfn通過單位球面表示法圖形化考察（超）極化率張量》（http://www.shanxitv.org/547）。
3 使用Multiwfn解析polar關鍵詞的輸出
對于上一節討論的各種情況，polar關鍵詞輸出的信息都能被Multiwfn所解析。注意route section一定要寫上#P，否則輸出文件將無法被正確解析。
啟動Multiwfn后，先輸入Gaussian的polar任務的輸出文件的路徑，比如C:\CH3NH2_polar.out，然后進入主功能200，選7，之后會看到下面的菜單：
-3 Set unit in the output, current: a.u.
-2 Set output destination, current: Output to screen
-1 Toggle loading frequency-dependent result for options 1 and 7, current: No
0 Return
1 "Polar" + analytic 3-order deriv. (HF,pure/hybrid DFT,semi-empirical)
2 "Polar" + analytic 2-order deriv. (MP2,double-hybrid DFT,TDDFT,CIS...)
3 "Polar=Cubic" + analytic 2-order deriv.
4 "Polar" + analytic 1-order deriv. (CISD,QCISD,CCSD,MP3,MP4(SDQ)...)
5 "Polar=DoubleNumer" or "Polar=EnOnly" + analytic 1-order deriv.
6 "Polar" + energy only (CCSD(T),QCISD(T),MP4(SDTQ),MP5...)
7 "Polar=gamma" + analytic 3-order deriv. (HF,DFT,semi-empirical)
選項1~7對應上一節提到的不同情況，當前你在什么情況下計算的就選哪種，然后Multiwfn就會輸出分析的結果。比如你用了CCSD結合polar=doublenumer就應當選5。之所以要對此進行選擇是因為不同情況下Gaussian的輸出格式相差極大，Multiwfn必須知道你在什么情況下計算的，才能相應地對輸出內容進行解析。選項1~7都可以給出α，只有1、3、5才會給出β，只有7才會給出γ。
如果在HF/DFT/半經驗計算時用了polar CPHF=rdfreq或polar=DCSHG做了含頻（超）極化率計算，或者用的是polar=gamma，默認情況下只會解析靜態的（超）極化率。如果想讓Multiwfn分析含頻的，應當先選-1將其狀態切換為Yes。當Gaussian輸入文件里定義了多個頻率時，如0.1,0.22,0.35，用戶可以選擇分析哪個頻率下的結果。對于β，用戶也可以選擇分析β(-ω;ω,0)還是β(-2ω;ω,ω)，注意只有寫了DCSHG時選擇β(-2ω;ω,ω)才有意義。對于γ，可以選擇分析γ(-2ω;ω,ω,0)還是γ(-ω;ω,0,0)。
默認情況下輸出的單位是a.u.，如果選-3，也可以切換為SI或esu單位。默認情況下輸出信息會顯示到屏幕上，也可以選-2來更改為輸出到當前目錄下的polar.txt當中。
如果你使用的是DFT或MP2，且對當前體系又要做振動分析又要算極化率，那么沒必要分別做freq和polar任務，只要做一次freq即可。用Multiwfn載入freq任務的輸出文件，在Multiwfn此文介紹的界面里選擇2 "Polar" + analytic 2-order deriv.即可提取極化率。
3.1 解析靜態極化率和第一超極化率一例
下面就以分析甲胺在#P PBE1PBE/aug-cc-pVTZ polar關鍵詞下得到的輸出文件為例進行說明。
當前例子用的是PBE1PBE，是DFT中常用的一種泛函，故對應的是前述第一種情況，因此當我們看到上面的菜單時應選1，然后馬上會看到α和β的分析結果，下面將對各項的含義進行說明。我們先看α部分
Dipole moment:
X,Y,Z= -0.506266 0.150290 0.000000 Norm= 0.528103
Static polarizability:
XX= 25.04620
XY= -0.08830
YY= 28.89110
XZ= 0.00000
YZ= 0.00000
ZZ= 24.43580
Isotropic average polarizability: 26.12437
Isotropic average polarizability volume: 3.871231 Angstrom^3
Polarizability anisotropy (definition 1): 4.18643
Polarizability anisotropy (definition 2): 4.18363
Eigenvalues of polarizability tensor: 24.4358 25.0442 28.8931
Polarizability anisotropy (definition 3): 4.15314
如前所述，默認情況下所有輸出的量都是以a.u.為單位。
極化率張量是個對稱的3*3矩陣，故所以只有6個獨立的矩陣元列了出來。每個體系的極化率的大小一般只用一個數值<α>來表示，也就是上面給出各項同性平均極化率，它是對XX、YY、ZZ取平均得到的。極化率也經常通過極化率體積來表示，上面也給出了。極化率并不是各項同性的，體系對不同方向射來的電場的響應是不同的。如上可見，Multiwfn也輸出了極化率的各向異性值Δα，數值越大各向異性越強，球對稱體系顯然其值為0。各向異性可以以三種方式定義，分別如下所示，詳見Multiwfn手冊3.200.7節
ε代表極化率張量的本征值，它也被Multiwfn輸出了出來。最常見的衡量極化率的各向異性程度是定義2，而定義3的形式一般是在衡量磁屏蔽張量各向異性時用的。
下面是第一超極化率部分的解析結果
Static first hyperpolarizability:
XXX= 13.173100
XXY= 4.937890
XYY= 13.647700
YYY= -35.484000
XXZ= 0.000000
XYZ= 0.000000
YYZ= 0.000000
XZZ= -2.763570
YZZ= 1.342190
ZZZ= 0.000000
Beta_X= 24.05723 Beta_Y= -29.20392 Beta_Z= 0.00000
Magnitude of first hyperpolarizability: 37.836745
Projection of beta on dipole moment: -31.373491
Beta || : -18.824094
Beta ||(z) : 0.000000
Beta _|_(z) : 0.000000
Gaussian的β的輸出有個眾所周知的問題，也就是符號是反的，即每個β元素的數值都應該乘以-1，至少直到G09 D.01版依然存在這個問題。有文章專門指出了這點，見《Gaussian程序計算的一階超極化率的符號問題》（四川師范大學學報（自然科學版）,33,228）。這個問題在Multiwfn中已經被考慮進去了，Multiwfn輸出的β的符號都是正確的。
β有3*3*3=27個元素。靜態β滿足Kleinman對稱性，三個標號可以隨意置換而完全不影響結果，比如xyy=yxy=yyx，只有10個元素是唯一的，所以Multiwfn也只輸出了10個元素。
Beta_X、Beta_Y、Beta_Z衡量的是β在X、Y、Z方向上的分量。Magnitude of first hyperpolarizability體現了β的整體的大小。其定義分別如下所示
只有β在偶極矩方向上的投影值才是可以通過電場誘導二次諧波產生實驗(EFISH)測量出來的，和實驗值對比通常比的都是這個值。Multiwfn給出了β在偶極矩方向上的投影值（β_prj，一些文中也寫為β_vec）。另一個文獻中常見的量是Beta ||，||代表平行于偶極矩，它正比于β_prj
文獻中常見的Beta ||(z)和Beta _|_(z)的定義如下所示，分別衡量的是平行和垂直于Z方向上的β的值
當偶極矩方向和Z方向一致時，Beta_Z=β_prj，Beta ||(z)=Beta ||。
值得一提的是，默認情況下，Gaussian會將輸入文件里體系的坐標調整到標準朝向下，這通常會令體系發生旋轉。Gaussian計算出的，也即Multiwfn所解析的，都是標準朝向下的結果，因此X、Y、Z可能和你的輸入文件里的笛卡爾軸方向不對應。如果想避免自動調整到標準朝向，用nosymm關鍵詞。
3.2 解析動態極化率和第一超極化率一例
此例使用Multiwfn解析polar關鍵詞產生的動態（超）極化率。輸入文件如下，分別計算外場為0.07和0.1 a.u.的情況。
#p PBE1PBE/aug-cc-pVTZ polar CPHF=RdFreq
test
0 1
C -0.55391731 0.43227932 0.14513431
H -1.29461381 -0.13844703 0.71478821
H -0.40071988 1.37810560 0.67504812
H -0.99858078 0.66864399 -0.83680924
N 0.70741247 -0.31205725 0.11165308
H 0.57479710 -1.18908235 -0.38659422
H 1.39864521 0.20863272 -0.42322026
0.07 0.1
令Multiwfn載入輸出文件，并進入主功能200里的子功能7后，先選-1切換成解析動態值，然后再選1。然后程序問載入哪個頻率下的結果，此例我們選2，即w=0.07的情況。然后我們得到了w=0.07下的極化率輸出。之后程序問分析β(-ω;ω,0)還是β(-2ω;ω,ω)的結果。由于沒寫DCSHG，我們只能選擇β(-ω;ω,0)。輸出如下
Frequency-dependent first hyperpolarizability Beta(-w;w,0)
XXX= 16.334900
XYX= YXX= 7.369050
YYX= 16.593300
XZX= ZXX= 0.000000
YZX= ZYX= 0.000000
ZZX= -2.918320
XXY= 8.947060
XYY= YXY= 15.473800
YYY= -40.341800
XZY= ZXY= 0.000000
YZY= ZYY= 0.000000
ZZY= 1.240500
XXZ= 0.000000
XYZ= YXZ= 0.000000
YYZ= 0.000000
XZZ= ZXZ= -2.796870
YZZ= ZYZ= 1.609830
ZZZ= 0.000000
Beta_X= 29.34451 Beta_Y= -30.96003 Beta_Z= 0.00000
Magnitude of first hyperpolarizability: 42.657048
Projection of beta on dipole moment: -36.941911
Beta || : -22.165146
Beta ||(z) : 0.000000
Beta _|_(z) : 0.000000
和上一例相比，可見在含頻外場下β值增大了。各項的含義和上例一致，不再累述。值得一提的是對于beta(-ω;ω,0)這種情況，只有i,j,k三個標號中的i,j可以互換，比如xyy=yxy≠yyx，張量的對稱性沒有靜態β的那么高，故有18個獨立的元素。不過，當外場頻率不是很大時，Kleinman對稱性是可以近似滿足的，從上面的結果可見，諸如xyx和xxy相差不算很大。
假設我們要看w=0.1時的β(-ω;ω,0)，那么再次進入此功能，然后選1-3-1即可。
基于當前輸出文件沒法查看β(-2ω;ω,ω)，因為當前例子沒寫DCSHG。要看它的話應把polar后面加上=DCSHG重算一次，然后再進入此功能，并選擇分析β(-2ω;ω,ω)情況的數據。
3.3 解析第二超極化率一例
這里我們對NH3計算第二超極化率γ。Gaussian輸入文件是Multiwfn程序包里的examples\polar\NH3_gamma.gjf，內容如下
#p PBE1PBE/daug-cc-pVTZ polar=gamma
[空行]
pbe1pbe/aug-cc-pVTZ opted
[空行]
0 1
N 0.00000000 0.00000000 0.11367600
H 0.00000000 0.93902800 -0.26524400
H -0.81322200 -0.46951400 -0.26524400
H 0.81322200 -0.46951400 -0.26524400
[空行]
532nm 680nm
想高精度計算小體系的γ，基組必須帶非常充足的彌散函數，因此此例用了每個角動量加兩層彌散的d-aug-cc-pVTZ，此基組非常昂貴。另外計算γ本身也很昂貴，所以哪怕是小體系，想高精度計算γ也需要給力的機子。polar=gamma關鍵詞要求既計算靜態的也計算動態的γ，如當前輸入文件所示，對于動態的γ我們要計算的是532 nm和680 nm外場下的情況。
此例我們要解析γ(-2w;w,w,0)在680 nm下的結果。啟動Multiwfn然后輸入
examples\polar\NH3_gamma.out
200
7 //解析Gaussian的(超)極化率的計算輸出
-1 //切換為載入動態(超)極化率（不選這個的話之后解析的將是靜態的γ）
7 //對應polar=gamma的情況
2 //對應680 nm的情況
2 //解析γ(-2w;w,w,0)形式的動態第二超極化率
輸出信息如下
XXXX= 3091.740000
YXXX= 0.000000
ZXXX= 0.000040
XYXX=XXYX= 0.000000
...略
YZYZ=YYZZ= 3677.510000
ZZYZ=ZYZZ= -0.446631
XZZZ= 0.000000
YZZZ= -0.400105
ZZZZ= 18072.400000
Magnitude of gamma: 5998.270291
X component of gamma: 1574.323333
Y component of gamma: 1574.284000
Z component of gamma: 5569.774000
Average of gamma (definition 1), gamma ||: 8718.3813
Average of gamma (definition 2): 8737.028000
gamma _|_: 2875.049333
上面每一項的定義筆者就不在這里介紹了，因為在Multiwfn手冊3.200.7節里已經寫得非常詳細了，這些量都是考察γ時經常要涉及的。
如屏幕上的提示所示，上面給出的gamma信息對應的是輸入朝向下的，而之前提到的Multiwfn給出的alpha和beta都是標準朝向下的，這個差異要注意。