使用Multiwfn基于完全態求和(SOS)方法計算極化率和超極化率

使用Multiwfn基于完全態求和(SOS)方法計算極化率和超極化率

文/Sobereva @北京科音

First release: 2014-Apr-27  Last update: 2020-Jul-5
    

1 前言

完全態求和(SOS)是計算極化率和超極化率的常見方法。Multiwfn支持通過SOS計算靜態/含頻的極化率α、第一超極化率β、第二超極化率γ乃至第三超極化率δ。其中需要用到的各個態的能量及偶極矩，以及各個態之間的躍遷偶極矩，都可以利用Multiwfn基于Gaussian或ORCA的CIS、TDHF、TDDFT計算結果來產生。本文就介紹具體怎么實現。

雖然Gaussian直接可以用TD(SOS)關鍵詞基于TDHF、TDDFT的結果做SOS計算，但是只能給出極化率，顯然太局限了，而且對常用的CIS也不支持SOS計算。

Multiwfn可以在其主頁http://www.shanxitv.org/multiwfn上免費下載，請用最新版本。Multiwfn是最強大的波函數分析程序，此文介紹的只不過是它的一個附屬功能而已。如果對Multiwfn不熟悉，建議閱讀《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。

另外，在《使用Multiwfn分析Gaussian的極化率、超極化率的輸出》（http://www.shanxitv.org/231）一文中，介紹了怎么使用Multiwfn分析Gaussian的polar關鍵詞產生的極化率和超極化率的輸出，其中也介紹了一些（超）極化率的基本概念，有興趣可以看看。polar關鍵詞所用的導數方法是另一種最重要的計算（超）極化率的方法。

Multiwfn的SOS功能還可以做對于考察第一超極化率本質非常重要的雙能級分析和筆者提出的三能級分析，在另一篇博文有專門說明：《使用Multiwfn對第一超極化率做雙能級和三能級模型分析》（http://www.shanxitv.org/512）。

2 SOS方法

SOS的具體原理這里不討論，主要就是介紹一下實際的計算公式。α、β、γ、δ的SOS計算公式在J. Chem. Phys., 99, 3738 (1993)里都能找到。注意很多很多文獻、書籍里給出的SOS公式都是錯的。計算極化率和第一超極化率的公式如下

A,B,C...這樣的標號用來表示方向{X,Y,Z}。ω是外場能量，為0時對應的是靜態的（超）極化率。加和是對所有激發態加和，Δ代表激發態相對于基態的激發能。P表示對方括號里的項進行各種可能的置換。比如對于β，P的方括號里有三項，因此有3!=6種排列，要對這6種情況的結果加和。μ_A_ij代表i、j兩個態的躍遷偶極矩的A方向的分量。i=j時對應的就是第i個態的偶極矩，故μ00就是基態的偶極矩。δij是kronecker符號，i=j時為1，否則為0。

類似地，第二、第三超極化率的公式為
 

可見，做SOS計算并不困難，公式都是現成的，只要提供各個激發態的激發能、偶極矩，以及各激發態之間的躍遷偶極矩就行了。通過ZINDO、CIS、TDHF、TDDFT等電子激發方法都能產生這些量。CIS(D)也可以，用于SOS時躍遷偶極矩依然是用CIS的，但激發能是通過二階微擾修正的，以此更好地考慮電子相關效應。

整個SOS計算過程分為兩部分(1)電子激發計算(2)依據SOS公式進行循環累加。對于(2)部分，α和β計算本身幾乎不花時間，γ稍微耗點時間。對于δ，當考慮的態數很多時，從上式可見需要對激發態做四重循環累加，而且包含3^5=243個分量（盡管有些分量是相同的可以避免重復計算），故計算δ還是頗費時的。對于一般我們感興趣的α、β和γ，SOS的整個計算花費主要都是用(1)上，特別是對于大體系、高質量基組的情況。注意(2)部分耗時本身和基函數數目無直接關系，只取決于考慮的態數，但(1)的耗時與基函數數目直接相關。

ZINDO是半經驗方法，計算很快，SOS/ZINDO組合很廉價，經常用在計算有機大共軛體系NLO性質上。而原理上更精準的CIS/TDHF/TDDFT顯然就耗時得多的多了。SOS原則上要對所有態進行求和，實際研究中雖然用不著考慮所有態，但一般也得在電子激發過程中計算40~120個態，遠高于一般研究電子激發問題要求的態數。求解的態數越多，CIS/TDHF/TDDFT就明顯越耗時。本來這類從頭算方式的電子激發計算用在很大體系就比較困難，再加上為了做SOS又要計算那么多態，而且較準確計算β，特別是γ，又需要彌散函數豐富的較大基組，SOS結合CIS/TDHF/TDDFT能應用的尺度是很有限的。比起SOS，如果所有的導數都能解析地計算的話，用導數方法計算（超）極化率是更好的選擇。不過支持高階解析導數從編程角度來講很困難，Gaussian能做到三階解析導數（對應β）的方法只有HF、DFT和半經驗，它們不支持四階解析導數來由此產生γ。雖然基于三階解析導數自行做一次有限差分也能得到靜態的γ，但是沒法得到含頻的γ，這就不得不用SOS了。另外，只要有了SOS計算所需的信息，每次SOS的計算很快，只不過是做簡單的循環及加減乘除而已，因此可以便利地研究α、β和γ隨頻率的變化，這是SOS的一個鮮明的優點。

3 Multiwfn的SOS計算功能

Multiwfn的SOS功能對應于主功能200里的子功能8，能夠計算極化率和第一、二、三超極化率，計算效率很高，被充分地并行化（機子是幾個核建議就在settings.ini文件里把并行線程數nthreads設成幾）。而且還能十分方便地計算α、β、γ隨著考慮的態數發生的變化以檢驗SOS的收斂性，以及α、β、γ隨著外場頻率的變化而發生的變化。

Multiwfn做SOS時要讀取的信息包括基態和激發態偶極矩、激發能，以及所有態間的躍遷偶極矩。這些信息可以從兩種渠道讀取
(1) 直接從Gaussian的ZINDO/CIS/TDHF/TDDFT輸出文件中讀取。不幸的是，TDHF/TDDFT任務的輸出文件不包含激發態之間的躍遷偶極矩，而CIS和ZINDO雖然有alltransitiondensities關鍵詞可以輸出激發態之間的躍遷偶極矩，但是并沒有辦法一次性給出每個激發態的偶極矩（除非改Gaussian代碼，比較復雜），而這些信息都是計算超極化率要用到的。故直接讀取Gaussian輸出文件，獲得的信息只夠計算極化率的，故此時Multiwfn也只輸出極化率。

(2)文本文件。用戶可以把各種各樣量化程序產生的SOS計算所需信息按照類似下面的格式寫到一個文本文件里，然后在Multiwfn啟動時讀取之并做極化率和超極化率計算。下面這個簡單例子假設僅考慮兩個激發態
2      // 激發態數
1  1.1        // 第1個激發態，其序號以及激發能(eV)
2  3.2
0 0  0.845 0.2 0.4    // 0代表基態，此項的含義是基態偶極矩的XYZ分量(a.u.)
0 1  0.231 0.3 0.7    // 基態到第1激發態的躍遷偶極矩
0 2  0.112 0.564 0.21
1 1  0.021 0.465 0.0    // 第1激發態的偶極矩
1 2  0.001 0.3 0.11     // 第1到第2激發態的躍遷偶極矩
2 2  0.432 0.14 0.42    // 第2激發態的偶極矩
可見，這種文本文件的意義很容易理解，也很容易編寫。即先給出激發態數，然后給出每個激發態的激發能，然后做i;j>=i循環，給出各個態的偶極矩以及態之間的躍遷偶極矩即可。格式是自由格式，也就是小數位數多少都可以。
（如果你只想計算極化率，那么1 1那行以及之后的行都可以不寫，因為這部分信息用不到，但此時需要把第一行寫為負數，例如上例寫為-2，這樣Multiwfn才能知道這一點）

Multiwfn可以直接基于Gaussian或ORCA的CIS/TDHF/TDDFT輸出文件計算出所有態之間的躍遷偶極矩以及每個態的偶極矩，這在《使用Multiwfn計算激發態間的躍遷偶極矩和各個激發態的偶極矩》（http://www.shanxitv.org/227）當中已經詳細介紹過，利用這些信息就可以用SOS方法計算極化率和超極化率了。為了方便起見，在Multiwfn中這些信息可以直接輸出成上述格式的文本文件，用于Multiwfn的SOS計算之目的。

下面就介紹一個完整的例子，利用Multiwfn基于Gaussian的CIS計算結果，通過SOS方法得到氨分子的極化率和第一、第二超極化率。對于SOS之目的，用TDHF相比CIS改進不大，而耗時增加了不少。TDDFT也同樣比CIS更耗時得多，用的雜化泛函的雜化成份越小，通常超極化率的計算結果比CIS大得越多（這是因為激發能也往往越低，使SOS公式的分母因此變小）。至于這令結果更準了還是更差了，不好說。

4 實例：NH3的SOS/CIS計算

本例計算使用Gaussian 09W B.01。

輸入文件的開頭部分如下，完整的輸入文件是Multiwfn的examples目錄下的NH3_SOS.gjf。
%chk=C:\gtest\NH3_SOS.chk
# CIS(nstates=150)/gen IOp(9/40=5)
[空行]
CCSD/def2-TZVPPD opted
[空行]
0 1
 N                  0.00000000    0.00000000    0.11461200
 H                  0.00000000    0.93646700   -0.26742900
 H                 -0.81100400   -0.46823400   -0.26742900
 H                  0.81100400   -0.46823400   -0.26742900
...[def2-TZVPPD基組定義部分略]

這里有幾個要點，首先是計算的態數。前面已經提到過，計算SOS需要算很多激發態。通常來說，越低的激發態對SOS的結果影響越大（但個別高階激發態可能也有很大影響），隨著考慮的態數增加結果會逐漸收斂。此例考慮150個態已經足夠了，結果肯定已經收斂了，下文還會專門檢驗這一點。做（超）極化率計算需要具有豐富彌散函數的基組才能得到較好定量結果，超極化率的階數越高，對彌散函數的數量和角動量要求越高。基組這里用的是def2-TZVPPD（JCP,133,134105），是在def2-TZVPP基組上為了得到較好極化率而增加了彌散函數的版本，計算第一超極化率也夠用了。由于這個基組不是Gaussian自帶的，所以從EMSL基組庫中拷出來，通過自定義基組方式使用。對于計算SOS之目的，寫上IOp(9/40=5)是必要的。因為默認情況下只有大于0.1的MO躍遷系數會被輸出到輸出文件中，而較小的都不輸出，這樣的話Multiwfn基于這些系數算出的各個態的偶極矩以及它們間的躍遷偶極矩結果將會不準確。IOp(9/40=x)的含義是將系數大于10^-x的組態都輸出出來，因此IOp(9/40=5)會把系數絕對值大于0.00001的組態都輸出，這種情況下Multiwfn給出的結果就足夠準確了（x設為大于5沒意義，因為Gaussian輸出時的格式只有5位小數）。計算過程中不要求必須用#P。

對于此例，假設要用TDDFT，就寫成比如# TD(nstates=150) CAM-B3LYP/gen IOp(9/40=5)即可。

將計算后得到的chk文件用formchk轉換成.fch文件，假設叫NH3_SOS.fch，并假設它和Gaussian輸出文件NH3_SOS.out都在C:\gtest目錄下。然后啟動Multiwfn，依次輸入
C:\gtest\NH3_SOS.fch
18    //電子激發分析功能
5     //計算態之間的躍遷偶極矩
C:\gtest\NH3_SOS.out
3     //將計算的結果，連同激發能等各種SOS所需信息，組合成Multiwfn的SOS功能所需的輸入文件
很快就計算完畢并在當前目錄下生成了SOS.txt文件，其格式和上一節介紹的一致。

重啟Multiwfn，依次輸入
SOS.txt
200    //主功能200（雜項功能的集合）
8      //用SOS方法計算（超）極化率
屏幕上輸出了每個態的激發能、基態偶極矩，并顯示出SOS計算菜單。如屏幕提示所述，一切輸出信息都以a.u.為單位。（超）極化率的a.u.、SI和esu單位的轉換關系可參見手冊3.200.7節末尾。

菜單中先選哪項都可以，沒有順序要求。假設我們先計算靜態極化率，選擇1，然后輸入0（即外場頻率為0），屏幕上立刻看到結果，如下所示。屏幕上的輸出都可以直接拷貝到文本文件里，如果不知道怎么做，參考手冊5.4節。

Multiwfn不僅輸出了極化率張量，還輸出了各種很重要的相關的量，它們的定義可參看前述的《使用Multiwfn分析Gaussian的極化率、超極化率的輸出》或手冊3.200.7節。從輸出中可見，各項同性平均極化率<α>的結果為14.59，和實驗值14.56(Mol. Phys.,33,1155)吻合得賊好。再來看看動態極化率。選1，輸入外場頻率0.0719，此時<α>的結果為14.86，比靜態的有所增大。輸入頻率的時候默認都是以a.u.為單位，但如果你想以nm為單位輸入頻率也可以，見屏幕上的提示。

接下來計算第一超極化率β(-(ω1+ω2);ω1,ω2)。我們主要感興趣的是β在偶極矩方向的分量β||，因為這是可以通過電場誘導二次諧波產生實驗(EFISH)測得的，因此可以用于檢驗結果的準確性。

先來計算β(0;0,0)。選2，然后輸入ω1和ω2的值，此時即輸入0,0。輸出如下

從輸出中看到β||值為-38.98。手頭沒有相應的實驗值，但可以和高精度CCSD(T)下計算的對比，其值為-34.3（JCP,98,3022），可見還是很接近的。

ω=0.0656下的β||(-2ω;ω,ω)是有實驗數據的，值為-48.9±1.2。我們也通過SOS方法算算看。選2，輸入0.0656,0.0656，從結果中看到β||值為-49.69，和實驗吻合得相當好。

不過需要提醒的是，盡管這一節的例子和實驗吻合得頗為不錯，但對于其它體系情況遠非總是這么樂觀。SOS/CIS對小分子往往明顯高估β，見A. Hernández-Laguna et al. (eds.), Quantum Systems in Chemistry and Physics, Vol. 1, p111）；也有明顯低估的時候，比如JCP,125,024101研究的兩種D-A型π共軛分子。

下面算一下第二極化率γ。選3，然后還是輸入外場頻率，這里輸入0,0,0以得到靜態γ。計算γ相對耗時一點，故程序默認不用所有的態，而允許你輸入考慮多少態。這里把所有態都考慮，故輸入150。程序輸出了γ張量的全部分量，以及一些相關信息。通常主要感興趣的是γ的平均值，結果為928.72。由于手頭沒有實驗值，精度如何不清楚，但可能不會太高，一方面是SOS/CIS在定量上本身算不上多精確，CIS對電子相關考慮不足，另一方面是基組用得不算大。要想很準確計算第二極化率，有時得用到諸如d-aug-cc-pVTZ這樣的程度（它比常用的aug-cc-pVTZ的彌散函數數目又多了一倍）。

雖然Multiwfn還能計算第三超極化率δ，但這里就不舉例了，因為計算時間會很長。而且δ并不重要，計算δ對方法、基組要求也都會極苛刻。

前面已經提到過，做SOS很重要的一個問題是所用的態數必須足夠多，SOS在原理上是要考慮所有激發態的。這也就是說，比如實際計算取了n個態，那么隨著考慮的態數的增多，考慮到n個態之前結果就必須已經基本收斂了。我們來作一下靜態極化率隨態數的變化圖。選擇5，然后輸入0，馬上在當前目錄下產生了alpha_n.txt文件。每一列的含義在屏幕上都已經清楚地提示了。將此文件直接拖入到Origin里，把第一列和第二列分別作為X和Y軸數據作曲線圖，就得到下面的左圖，展現了<α>隨著態數增加的收斂性。類似地，若要研究靜態β||的收斂性，選6，輸入0,0，對輸出的beta_n.txt里的第一列和第七列作成曲線圖，就得到了下圖的右圖。另外還可以用選項7研究γ的收斂性。

從上圖可見，無論是<α>還是β||，在100個態的時候都已經收斂了，繼續增加態數結果變化很小，因此此例用了150個態是足夠多的，或者說只取150個態不會造成明顯的態截斷誤差。如果只是要得到定性正確的結果，實際上取70個態也勉強夠了。SOS的收斂性和體系特征關系密切，需要用多少得實際進行檢驗，很多大體系往往取40、50個態就足夠了。

用Multiwfn還可以方便地研究（超）極化率隨外場頻率的變化。如果要研究極化率，選15，然后輸入頻率的起始值、終止值和步長就行了，每個頻率下的極化率都會被計算并輸出到當前目錄下的alpha_w.txt里，然后也是可以直接用origin對相應的列進行作圖。對于研究β和γ隨頻率的變化，可以分別用選項16和17實現。由于β和γ分別涉及到兩個和三個頻率，存在倍頻、差頻、和頻等情況，比較復雜，用戶需要自行編寫一個文本文件讓Multiwfn載入，每一行對應一對兒要計算的外場頻率。比如我們這里研究β(ω;ω,0)隨著ω從0增加到0.2過程中的變化，步長為0.01。為了方便，我們用excel，在A1~A21單元格通過拖拉產生0.0,0.01,0.02...0.2序列，并令B1~B21等于0。然后選擇“另存為”，格式選擇文本文件（制表符分隔），保存到比如C:\freqlist.txt里面，文件內容應當是這樣
0 0
0.01 0
0.02 0
...
0.19 0
0.2 0
在Multiwfn里選16，然后輸入C:\freqlist.txt。很快就算完了，結果輸出到了當前目錄下beta_w.txt下面。每一列含義在屏幕上也明示了。我們將此文件拖進Origin，對第一列(ω)和第8列(β||)進行作圖，結果如下

程序還同時在當前目錄下輸出了beta_w_comp.txt，里面是beta的所有分量隨頻率的變化。

Multiwfn的SOS功能還可以對β(-(ω1+ω2);ω1,ω2)的兩個入射光頻率做二維掃描，得到下面這種圖，對于討論一些問題很有益，詳見此帖http://bbs.keinsci.com/thread-15455-1-1.html。

Multiwfn的SOS功能的界面里還有很多其它選項，限于篇幅就不一一介紹了，一看相應的選項文字和屏幕上的提示就秒懂，請結合手冊里的說明自行嘗試。