使用Multiwfn計算激發態之間的密度差

使用Multiwfn計算激發態之間的密度差

文/Sobereva @北京科音

First release: 2018-Jul-28   Last update: 2022-Jul-14

1 前言

使用Multiwfn計算密度差、繪制各種密度差圖形非常簡單，在《使用Multiwfn作電子密度差圖》（http://www.shanxitv.org/113）中有充分示例，其中也包括繪制激發態與基態間的密度差，這對于討論電子躍遷過程中的電子轉移特征很有意義。對Multiwfn稍有了解的人肯定也都明白，通過Multiwfn同樣也可以繪制激發態之間的密度差，只要讓Gaussian等量化程序輸出記錄了相應兩個激發態的波函數的文件，然后按照常規步驟繪制密度差即可。

但是，有時由于特殊需要，我們要對一大批激發態兩兩之間計算密度差，如果用量子化學程序對每個激發態都做一次計算來產生記錄了其波函數的文件，顯然是很繁瑣的事情；即便可以用腳本自動化實現，但總耗時也比只做一次電子激發計算高很多。好在Multiwfn可以讀取Gaussian等程序的激發態計算任務的輸出文件，一次性就能把各個激發態的.mwfn文件全都產生出來。.mwfn是Multiwfn支持的記錄波函數信息的一種格式，見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379），對于TDDFT等方法計算的激發態，這樣得到的.mwfn文件里記錄的是激發態的自然軌道。之后想對哪兩個激發態之間繪制密度差就利用哪兩個態的.mwfn文件即可，非常方便，本文就結合Gaussian示例一下操作。更多關于此功能的說明可參見Multiwfn手冊3.21.13節。

另外，有了激發態的.mwfn文件后，將之載入Multiwfn后還可以對激發態做各種波函數分析，比如計算電荷分布、鍵級等等，分析實例可參看《在Multiwfn中基于fch產生自然軌道的方法與激發態波函數、自旋自然軌道分析實例》（http://www.shanxitv.org/403）等文章。

Multiwfn可以從其主頁http://www.shanxitv.org/multiwfn上免費下載，不了解Multiwfn的讀者建議參看《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）。Gaussian用的是G16W A.03。本文涉及的文件可在此下載：file.zip

2 實例：繪制N-苯基吡咯的激發態之間的密度差

首先使用以下輸入文件，在PBE0/6-31G*級別下做TDDFT計算，得到最低五個單重態激發態的信息。IOp(9/40=4)是為了讓Gaussian輸出所有絕對值大于1E-4的組態系數，這些組態系數連同基態分子軌道將被用于構建激發態自然軌道。不寫IOp(9/40=4)的話，默認情況下只輸出絕對值大于0.1的組態系數，此時Multiwfn產生的激發態自然軌道將很不準確。

%chk=C:\N-phenylpyrrole.chk
# PBE1PBE/6-31g(d) TD(nstates=5) IOp(9/40=4)

B3LYP/6-31G* opted

0 1
 C                 -0.00000000    1.12162908    1.82507914
 C                 -0.00000000    0.71310006    3.13457424
 C                 -0.00000000   -0.71310006    3.13457424
 C                 -0.00000000   -1.12162908    1.82507914
 N                 -0.00000000    0.00000000    1.00627307
 H                 -0.00000000    2.11821416    1.41565311
 H                 -0.00000000    1.36582111    3.99707730
 H                 -0.00000000   -1.36582111    3.99707730
 H                 -0.00000000   -2.11821416    1.41565311
 C                  0.00000000   -0.00000000   -0.41287403
 C                  0.00000000   -1.20626109   -1.13056908
 C                  0.00000000    1.20626109   -1.13056908
 C                  0.00000000   -1.20041809   -2.52387419
 H                  0.00000000   -2.15723816   -0.61186005
 C                  0.00000000    1.20041809   -2.52387419
 H                  0.00000000    2.15723816   -0.61186005
 C                  0.00000000   -0.00000000   -3.23381724
 H                  0.00000000   -2.14896316   -3.05414223
 H                  0.00000000    2.14896316   -3.05414223
 H                  0.00000000   -0.00000000   -4.31975533

將得到的chk用formchk轉換為fch文件。然后啟動Multiwfn，依次輸入

N-phenylpyrrole.fch  //剛得到的fch文件。里面記錄的是基態的DFT軌道
18   //電子激發分析
13   //產生激發態的.mwfn文件
N-phenylpyrrole.out  //剛得到的輸出文件
1-3   //假設當前我們只需要產生前三個激發態的.mwfn文件
馬上，當前目錄下就產生了NO_0001.mwfn、NO_0002.mwfn、NO_0003.mwfn，分別是記錄了第1、2、3激發態自然軌道的.mwfn文件。

假設我們目前想考察第3激發態（S3）與第1激發態（S1）之間的密度差的等值面圖，就重啟Multiwfn，然后輸入
NO_0003.mwfn
5  //計算格點數據
0  //自定義操作
1  //之后有一個文件會作用到一開始載入的文件上
-,NO_0001.mwfn   //令NO_0003.mwfn的屬性減去NO_0001.mwfn的屬性
1   //屬性設為電子密度
2   //中等質量格點
-1   //觀看等值面圖
在圖形窗口中把等值面數值設小到0.005，看到下圖。綠色區域為正，藍色區域為負，分別代表S3比S1密度多和少的區域。

實際上之后還可以做更多的分析，比如可以按照手冊4.13.6節示例的將密度差轉化為電荷位移曲線、按照手冊4.18.3節示例的將密度差的正、負部分平滑化為高斯分布并計算電荷轉移長度等等，這里就不再多述了。

Multiwfn如上方式產生的激發態的.mwfn文件對應的是非弛豫的激發態密度，這沒有弛豫的激發態密度那么真實，但好處是不需要額外耗時。如果非要用弛豫的激發態密度來求差，那就不得不對每個激發態分別做計算來產生波函數文件了，而且每次耗時和計算激發態能量梯度的耗時相仿佛，故代價較高。使用比如# PBE1PBE/6-31g(d) TD(nstates=5,root=3) out=wfn關鍵詞所產生的.wfn文件記錄的就是第3激發態的弛豫密度對應的自然軌道。關于弛豫和非弛豫密度的更多信息見http://bbs.keinsci.com/thread-5738-1-1.html。