在Multiwfn中基于fch產生自然軌道的方法與激發態波函數、自旋自然軌道分析實例
在Multiwfn中基于fch產生自然軌道的方法與激發態波函數、自旋自然軌道分析實例
文/Sobereva @北京科音
First release: 2018-Jan-8 Last update: 2020-Nov-6
0 前言
在《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）和Multiwfn手冊第四章開頭都說過，對于Gaussian用戶，如果要在Multiwfn中分析后HF波函數、TDDFT激發態波函數或者雙雜化泛函的波函數，需要產生相應級別的自然軌道并讓Multiwfn載入。自然軌道可以讓Gaussian以.wfn/wfx形式輸出，也可以讓Gaussian將之存到.fch里。前者很容易，但此類輸入文件在Multiwfn里沒法進行需要利用基函數信息的分析（如Mayer鍵級、Mulliken電荷）；后者雖然沒這個限制，但是導出過程很繁瑣，需要Gaussian中運行兩步，要用的關鍵詞容易忘，而且之后還得在.fch開頭自行填上saveNO。fch文件里是包含密度矩陣信息的，自然軌道也正是通過對角化密度矩陣得到的。但fch里的密度矩陣不會被Multiwfn所直接利用，Multiwfn從輸入文件里只會讀取軌道信息。從2018-Jan-7更新的Multiwfn開始，為了方便用戶，筆者在Multiwfn中添加了一個功能，可以讀取fch里的密度矩陣并產生各類自然軌道。本文目的是說一下相關細節，演示一下怎么用，并順帶簡要講一下激發態波函數的分析和自旋自然軌道的分析。
1 關于fch里的密度矩陣和自然軌道
在chk里，以及轉換出的fch文件里，至少包含了SCF級別的密度矩陣，對于后HF、雙雜化泛函、TDDFT等任務，還會包含其它類型密度矩陣（前提是用了density或某些特殊關鍵詞，而且當前理論方法在Gaussian中有解析梯度）。在fch文件里搜索Density就能知道包含哪些。這里羅列一下常見情況：HF、DFT、CASSCF任務：SCF
CIS、TDDFT任務：CI、CI Rho(1)
CCSD任務：CC
MP2、雙雜化泛函任務：MP2
MP3任務：MP3
MP4SDQ任務：MP4
如果你的體系是開殼層（對于CIS/TD來說指的是參考態），還會包含自旋密度矩陣。比如，對一個三重態體系，如果你用的是比如# MP2/cc-pVDZ density關鍵詞，則fch文件里會有這幾個字段：
Total SCF Density：HF級別的總密度矩陣
Spin SCF Density：HF級別的自旋密度矩陣
Total MP2 Density：MP2級別的總密度矩陣
Spin MP2 Density：MP2級別的自旋密度矩陣
如果體系是閉殼層，就沒有Spin SCF Density和Spin MP2 Density了。
所謂總密度矩陣，就是指Alpha和Beta密度矩陣的加和；而自旋密度矩陣，就是Alpha密度矩陣減去Beta密度矩陣。
Multiwfn能產生的自然軌道有三類，和密度矩陣的關系如下：
(a)自然軌道(Natural orbital, NO)：對總密度矩陣對角化得到。本征值是占據數，原則上范圍是0.0~2.0，本征矢就是自然軌道系數。這種自然軌道也被稱為Spatial natural orbital，因為不含自旋信息。如果體系是開殼層的，這種自然軌道也被確切叫做unrestricted natural orbital (UNO)。
(b)Alpha自然軌道、Beta自然軌道：分別對Alpha密度矩陣和Beta密度矩陣對角化得到，分別描述Alpha和Beta電子，原則上占據數在0.0~1.0之間。這類自然軌道也被稱為自然自旋軌道(natural spin orbital, NSO)。
(c)自旋自然軌道(Spin natural orbital, SNO)：對自旋密度矩陣對角化得到，占據數在-1.0至1.0之間。占據數為正和為負的SNO分別描述Alpha單電子和Beta單電子。
2 實例1：分析甲醛TDDFT計算的S1態波函數
運行以下輸入文件，對甲醛在TD-PBE0/6-31G*級別下做電子激發計算，并且把S1態密度矩陣存到chk里（因為TD里的root選項默認為1）。%chk=C:\H2CO.chk
# PBE1PBE/6-31G* TD density
B3LYP/6-31G* opted
0 1
C 0.00000000 0.00000000 -0.56221066
H 0.00000000 -0.92444767 -1.10110537
H -0.00000000 0.92444767 -1.10110537
O 0.00000000 0.00000000 0.69618930
用formchk將chk轉換為fch，里面此時會有Total SCF Density、Total CI Rho(1) Density、Total CI Density三個字段，分別對應于基態DFT密度矩陣、TDDFT非弛豫的S1態密度矩陣（可無視）、TDDFT弛豫的S1態密度矩陣。
下面我們來考察一下S1態的C-O Mayer鍵級以及O的ADCH原子電荷。啟動Multiwfn，依次輸入
H2CO.fch
200
16 //基于fch里的密度矩陣產生自然軌道
CI //載入的密度矩陣級別
由于輸入的是CI，所以Multiwfn從fch中讀取了Total CI Density字段記錄的S1態總密度矩陣，并立刻產生了相應的S1態自然軌道，占據數直接輸出在屏幕上了（占據數稍微超過理論上0.0~2.0范圍是正常情況，不用糾結）：
Occupation numbers:
2.033734 2.002238 2.000812 2.000459 2.000148 2.000000
2.000000 0.996739 0.970938 0.000094 0.000062 0.000049
0.000019 0.000005 0.000002 0.000000 0.000000 0.000000
0.000000 0.000000 -0.000000 -0.000000 -0.000000 -0.000000
-0.000000 -0.000000 -0.000000 -0.000000 -0.000238 -0.000366
-0.000454 -0.000828 -0.001168 -0.002245
此時，內存里的基函數信息已經被更新成了對應S1自然軌道的情況了，做各種基于基函數的分析的結果將是對應S1態的了。我們還要計算ADCH電荷，這是要用GTF信息的，但目前GTF系數還是對應基態DFT軌道的的（因為一開始載入的.fch里的軌道是基態DFT的），且現在軌道占據數已經變成了自然軌道占據數，所以之后做基于GTF信息的分析將沒有任何物理意義。為了能計算S1態的ADCH電荷，我們此時應當選y，讓程序把當前的波函數信息導出為當前目錄下的new.mwfn，并自動載入之，此時內存里的GTF信息也變成對應S1態的了。這個new.mwfn文件我們以后也再可以反復利用。
接著輸入以下內容
0 //退回主菜單
9 //鍵級計算
1 //Mayer鍵級
得知C-O鍵級為1.18。然后輸入
n //不導出鍵級矩陣
0 //退回主菜單
7 //布居分析
11 //ADCH電荷
1 //用內置的球對稱化的自由原子密度
得知O的ADCH電荷為-0.128。
我們可以對比一下基態的情況。重啟Multiwfn，載入H2CO.fch，此時內存里的基函數信息和GTF信息都是對應基態的，重做上述計算，發現C-O Mayer鍵級為1.98，O的ADCH電荷為-0.296。
之所以激發態的C-O鍵級比基態的小得多，激發態的O的電荷也比基態的正得多，是因為此激發是n->pi*激發，一方面削弱了C-O鍵，一方面一部分電子從氧轉移給了C。
3 實例2：分析丁烷雙自由基的自旋自然軌道
丁烷雙自由基體系在《CASSCF計算雙自由基以及雙自由基特征的計算》（http://www.shanxitv.org/264）中做了很詳細的討論，建議看看。本節的例子探究一下此體系的SNO軌道。
運行以下輸入文件
%chk=C:\C4H8.chk
#p UB3LYP/def2SVP guess=mix nosymm
ub3lyp/6-31g(d) opted
0 1
C -0.74400100 1.78566400 0.00000000
H -0.60282700 2.33865300 0.92499500
H -0.60282700 2.33865300 -0.92499500
C -0.74400100 0.30988100 0.00000000
H -1.25452600 -0.08746700 0.88463900
H -1.25452600 -0.08746700 -0.88463900
C 0.74400100 -0.30988100 0.00000000
H 1.25452600 0.08746700 -0.88463900
H 1.25452600 0.08746700 0.88463900
C 0.74400100 -1.78566400 0.00000000
H 0.60282700 -2.33865300 -0.92499500
H 0.60282700 -2.33865300 0.92499500
將chk轉換成fch后，啟動Multiwfn，依次輸入
C4H8.fch
200
16
SCF //考察DFT級別的密度矩陣
此時程序問你，是產生空間自旋軌道（此時叫做UNO），Alpha和Beta各自的自然自旋軌道(NSO)，還是自旋自然軌道(SNO)。我們選3產生SNO，看到的占據數如下
0.957908 0.045983 0.043627 0.037531 0.021487 0.021280
0.020200 0.017324 0.008494 0.007915 0.006743 0.006539
0.000626 0.000604 0.000082 0.000082 0.000000 0.000000
...略
-0.000000 -0.000000 -0.000082 -0.000082 -0.000604 -0.000626
-0.006539 -0.006743 -0.007915 -0.008494 -0.017324 -0.020200
-0.021280 -0.021487 -0.037531 -0.043627 -0.045983 -0.957908
可見占據數為正的SNO當中只有一個數值較大的，即第一個（0.9579）；而占據數為負的SNO當中也只有一個數值較大的，即最后一個（-0.9579）。而且，我們也知道，由于此體系是個很理想的雙自由基體系，本來就應該有一個Alpha和一個Beta單電子。因此，未成對的Alpha電子和Beta電子分布特征通過考察這兩個軌道就足夠得知了。
我們選y產生new.mwfn并自動載入之，然后進入主功能0觀看這兩個軌道的特征。值得一提的是，new.mwfn被Multiwfn視為是開殼層的，因此形式上會有Alpha和Beta兩套軌道，而我們產生的SNO軌道只有一套，此時SNO軌道占用的是用于記錄Alpha軌道那部分空間（絕對不是說SNO的物理意義是Alpha軌道），而Beta軌道那部分沒有意義，占據數也都是0。我們在主功能0的界面里可以選左上角的Orbital info.，從文本窗口里找出占據數為0.957908和-0.957908的軌道，可知序號分別是1和96（其實96就是體系的基函數數目，因此都沒必要去列表里找，直接選就行），這兩個軌道的圖形如下：
可見，Alpha單電子主要分布在C1上，在C4-C7之間也有少量分布。Beta單電子主要分布在C10上，同樣在C4-C7之間也有少量分布。
如果我們想定量知道比如Alpha單電子在各個原子上的量，我們可以對主要表現Alpha單電子的SNO1做軌道成分分析。我們回到主菜單，依次輸入
8 //軌道成分分析
8 //Hirshfeld方式做軌道成份分析
1 //用內置的球對稱化的自由原子密度
1 //1號軌道
結果為
Atom 1(C ) : 70.148760%
Atom 2(H ) : 6.317230%
Atom 3(H ) : 6.317230%
Atom 4(C ) : 7.588237%
Atom 5(H ) : 1.254212%
Atom 6(H ) : 1.254212%
Atom 7(C ) : 5.065489%
Atom 8(H ) : 0.591654%
Atom 9(H ) : 0.591654%
Atom 10(C ) : 0.758644%
Atom 11(H ) : 0.056339%
Atom 12(H ) : 0.056339%
數值和我們從軌道圖形上看到的很好相符。
眾所周知，研究單電子分布最常用的方法是分析自旋密度。我們重新載入C4H8.fch，繪制一下自旋密度圖，步驟見《談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算》（http://www.shanxitv.org/353），圖像如下：
我們再用主功能15，在Hirshfeld劃分下計算自旋密度在各個原子空間中的積分值（原子自旋布居）：
Atomic space Value
1(C ) 0.70928760
2(H ) 0.04488357
3(H ) 0.04488357
4(C ) -0.01371421
5(H ) 0.00849956
6(H ) 0.00849956
7(C ) 0.01371421
8(H ) -0.00849956
9(H ) -0.00849956
10(C ) -0.70928760
11(H ) -0.04488357
12(H ) -0.04488357
無論是從等值面圖上，還是原子自旋布居上，都看到自旋密度分布特征基本等于主要描述Alpha單電子的SNO1和主要描述Beta單電子的SNO96特征的疊加。對于兩個SNO基本沒有重疊的體系兩端的亞甲基的碳，其自旋布居和它對SNO的貢獻量十分接近，而在C4、C7部分，由于兩個SNO重疊厲害，導致顯著的相互抵消，所以它們雖然對SNO貢獻不很小，但是自旋布居卻非常小。
可見，分析SNO的好處是可以分別考察未成對的Alpha和Beta電子原本是怎么分布的，且可以從軌道分析的角度研究；而考察自旋密度，則便于了解在Alpha和Beta單電子分布發生一定程度抵消后的單電子“凈”分布情況，分析也更為省事，不用討論多個軌道。兩種分析是互補的。