使用Multiwfn做自然躍遷軌道(NTO)分析

2022-Mar-1補充：Multiwfn已經支持結合CP2K程序對周期性體系做NTO分析，見《使用CP2K結合Multiwfn對周期性體系模擬UV-Vis光譜和考察電子激發態》（http://www.shanxitv.org/634）。

使用Multiwfn做自然躍遷軌道(NTO)分析

文/Sobereva@北京科音
First release: 2017-May-26  Last update: 2023-Feb-8

1 前言

以前筆者專門寫過一篇文章介紹NTO的原理和用處，見《躍遷密度分析方法-自然躍遷軌道(NTO)簡介》（http://www.shanxitv.org/91）。NTO在電子激發分析方面比較流行，已被Gaussian、ORCA、Q-Chem等諸多量子化學程序直接支持。Multiwfn程序中也支持NTO分析，相對于用那些量子化學程序自帶的NTO功能的好處在于在Multiwfn中做NTO分析超級方便，想分析哪個激發態就輸入哪個態的編號即可，輸出信息也特別易于理解。Multiwfn也能很便利地可視化NTO、分析NTO的軌道成份，甚至還可以計算NTO的軌道能量。而且利用筆者現成的shell腳本還可以實現一行命令就讓Multiwfn把所有激發態的NTO軌道全都產生出來，超級快捷省事。反觀用量子化學程序自帶的NTO，對多個激發態一一做NTO分析那就麻煩太多了，耗時也比用Multiwfn高得多得多，還很不靈活。所以對于做NTO分析的人，筆者強烈建議通過Multiwfn實現。

本文目的在于簡單介紹怎么在Multiwfn里做NTO分析，并示例如何快速批量做NTO分析。本文對應官網上最新版本的Multiwfn的情況。Gaussian用的是G09 D.01。Multiwfn可在其主頁http://www.shanxitv.org/multiwfn免費下載。對Multiwfn不了解的話建議參看《Multiwfn入門tips》（http://www.shanxitv.org/167）、《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）。NTO的原理在前述文章中已經介紹過，一些詳細細節在手冊3.21.6節也有介紹，因此在此文不再累述。另外，Multiwfn還支持一大堆其它非常重要的電子激發分析方法，總結見《Multiwfn支持的電子激發分析方法一覽》（http://www.shanxitv.org/437），強烈建議一看！

2 輸入文件

Multiwfn做NTO分析需要兩類文件。下面說TDDFT的時候一律也包括了CIS、TDHF、TDA-DFT的情況。
(1)含有基函數信息的文件，比如.fch、.molden、.gms等，詳見Multiwfn手冊2.6節，其中應記錄了TDDFT計算的參考態軌道。
(2)記錄TDDFT計算產生的含有電子激發信息的文件，因為Multiwfn做NTO要從中讀取激發態的組態系數。可以直接用Gaussian的CIS/TD/TDA關鍵詞的輸出文件（必須是對單個結構計算，不能是幾何優化等），也可以用滿足格式要求的文本文件，看手冊3.21.A節對格式的描述，這使得其它量化程序利用Multiwfn做NTO分析成為可能（自己寫個程序轉換一下激發態輸出信息格式即可）。

對于Gaussian用戶，簡而言之，就是用CIS、TD或TDA關鍵詞對某個結構做電子激發計算，將輸出文件以及同時產生的fch文件載入到Multiwfn里即可。注意計算時必須加上IOp(9/40=4)關鍵詞，這使得組態系數絕對值大于>0.0001的都能夠被輸出，而Gaussian默認只輸出組態系數絕對值大于>0.1的組態系數，會導致NTO分析的精度明顯不夠。用IOp(9/40=5)會使NTO分析更精確一丁點，但對于大體系，可能輸出文件會變得很大，NTO分析耗時也會增加，故沒必要。

Multiwfn也直接支持基于ORCA做NTO分析，用其輸出的.molden文件代替.fch文件，用ORCA的激發態任務輸出文件代替Gaussian的輸出文件即可。Multiwfn的NTO功能對輸入文件的要求詳見Multiwfn手冊3.21.A節。

Multiwfn還支持NTO軌道能量的計算，這還需要用戶提供含有Fock或Kohn-Sham矩陣的文件，格式說明見Multiwfn手冊3.300.6節，計算NTO軌道能量的完整例子見4.300.6節。

3 NTO分析實例：尿嘧啶

這里以尿嘧啶(uracil)作為例子演示用Multiwfn基于Gaussian的輸出做NTO分析。如果對Gaussian做電子激發計算的基本常識都沒有的話應當先看《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）。

下面是尿嘧啶的激發態計算的輸入文件，結構已經在B3LYP/6-31G*下對基態優化好。這個任務會在TD-PBE0/6-31G*級別下計算最低三個單重態激發態。
%chk=C:\uracil.chk
#p PBE1PBE/6-31G* TD IOp(9/40=4)

B3LYP/6-31G* opted

0 1
 C                 -1.23749815    0.36573215   -0.00000000
 C                  0.05337385   -1.70928465    0.00000000
 C                  1.24151979   -1.06801061    0.00000000
 C                  1.28200385    0.39168010    0.00000000
 H                 -0.02061908    2.00141365   -0.00000000
 H                 -2.02597565   -1.51440798   -0.00000000
 H                 -0.03161872   -2.79076675    0.00000000
 H                  2.18093078   -1.60364321    0.00000000
 N                 -0.00000000    0.98771295   -0.00000000
 N                 -1.14105931   -1.02649028   -0.00000000
 O                  2.28466327    1.08546823    0.00000000
 O                 -2.30362555    0.95179972   -0.00000000

計算完成后，把.chk轉換為.fch。此任務的輸出文件和相應的.fch文件已經在Multiwfn壓縮包自帶的examples\excit\NTO目錄下提供了。

我們將要對這個體系的S0->S3激發做NTO分析，在做分析之前，我們先看看不寫IOp(9/40=4)時候做這個TDDFT計算時輸出的S3的組態系數：
26 -> 30  0.54135
26 -> 31 -0.20634
28 -> 30 -0.15424
28 -> 31  0.36715
可見，沒有哪個MO對對激發起到主導作用，貢獻最大的也僅僅是0.54135^2*2*100%=58.5%。如果不會算貢獻值的話看《電子激發任務中軌道躍遷貢獻的計算》（http://www.shanxitv.org/230）。因此，對這個激發，顯然不能光靠只分析一對MO軌道來討論其特征，而同時也討論其它的MO對則會很麻煩，不好交代清楚。NTO分析對于解決這種困難很有效，因為把MO變換成NTO之后，往往就只有一對NTO對電子激發有很大貢獻了，屆時就只討論這一對NTO的特征就夠了。

下面我們開始對尿嘧啶的S0->S3做NTO分析。啟動Multiwfn，然后依次輸入
examples\excit\NTO\uracil.fch
18   //電子激發分析模塊
6   //NTO分析
examples\excit\NTO\uracil.out
3   //對第3個激發態作分析
馬上看到如下信息
Multiplicity of this excited state is   1
Excitation energy is   6.0180000 eV
There are     846 orbital pairs in this transition mode
The sum of square of excitation coefficients:  0.500203
The negative of the sum of square of de-excitation coefficients: -0.000196
The sum of above two values  0.500007
以上信息中顯示，這個激發態自旋多重度是1，激發能是6.018eV。Multiwfn共從輸出文件中讀入了846個組態系數，貢獻和為0.500007，和理想值0.5相差甚微，所以當前的NTO分析精度是足夠高的。

再往下會看到NTO對的本征值，為了簡明只列出了最大的10個本征值，從大到小排序。
The highest 10 eigenvalues of NTO pairs:
   0.865529    0.134025    0.000582    0.000121    0.000063
   0.000024    0.000016    0.000015    0.000007    0.000006
Sum of all eigenvalues:  1.000387

從以上數據，我們就能知道對當前S0->S3躍遷，貢獻最大的一個NTO對的貢獻值是86.55%。由于此數值已經很大，足夠認為產生了主導效應，因此我們只要分析這個NTO對所對應的兩個NTO的特征，就可以搞清楚S0->S3躍遷的主要特點了。

現在，Multiwfn會讓你選擇是否導出含有剛產生的NTO軌道的.fch或.molden或.mwfn文件，導出哪種格式都可以，之后都可以被Multiwfn載入用于觀看和分析軌道。這里我們選擇把NTO軌道導出成.fch文件，然后輸入導出的文件的路徑，比如輸入C:\nicomaki.fch。之后Multiwfn會重新載入examples\NTO\uracil.fch恢復一開始的狀態。之后如果比如你再想對S0->S2做NTO分析，那么可以再次進入NTO分析功能，然后輸入2。可見非常方便，輸入文件就是一套，想對哪個激發態分析NTO就選幾即可。

下面我們看看我們剛才對S0->S3產生的那些NTO。重啟Multiwfn，載入剛生成的C:\nicomaki.fch，之后進入主功能0查看軌道。如果不熟悉Multiwfn看軌道的功能可以參看《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）。我們在圖形界面的左上角選orbital info. - Show up to LUMO+10，此時從1號到LUMO+10軌道的信息就都輸出在文本界面了，軌道能量那一列此時已經不是軌道能量了，而是NTO本征值（NTO軌道沒有能量的概念）。這里把其中一部分信息截出來：
Orb:    27 Ene(au/eV):     0.000582       0.0158 Occ: 2.000000 Type: A+B
Orb:    28 Ene(au/eV):     0.134025       3.6470 Occ: 2.000000 Type: A+B
Orb:    29 Ene(au/eV):     0.865529      23.5522 Occ: 2.000000 Type: A+B
Orb:    30 Ene(au/eV):     0.865529      23.5522 Occ: 0.000000 Type: A+B
Orb:    31 Ene(au/eV):     0.134025       3.6470 Occ: 0.000000 Type: A+B
Orb:    32 Ene(au/eV):     0.000582       0.0158 Occ: 0.000000 Type: A+B
可見，占據的NTO29和非占據的NTO30組成了一個NTO對，從a.u.為單位的能量那一列可以看到，它們就是本征值為0.8655，即對S0->S3躍遷貢獻高達86.55%的那一對NTO。我們把這兩個軌道的圖在Multiwfn里顯示出來，如下所示

可以清楚地看出，NTO29可以認為是O12的孤對電子軌道，NTO30是六元環的pi*軌道，因此這個激發模式可以指認為n->pi*。

從NTO本征值上也可以看到貢獻次高的是NTO28-NTO31這一對，貢獻為13.4%，有興趣的話也可以繪制軌道圖形分析一下它們對應的特征。

很值得一提的是，Multiwfn有大量分析軌道的功能，對于任何軌道類型都可以用，也包括NTO。比如，我們要通過SCPA方法分析NTO29中O12的貢獻，那么在載入nicomaki.fch后，可以進入主功能8，選擇3，然后輸入29，立刻各個原子對此軌道的貢獻就知道了，其中O12貢獻為60.47%。關于軌道成份計算詳見《談談軌道成份的計算方法》（http://www.shanxitv.org/131）。再比如，我們還可以討論NTO29與NTO30之間的重疊程度以及質心距離，從而定量討論S0->S3激發的電荷轉移(CT)的程度，詳見《使用Multiwfn考察軌道間重疊程度和質心距離》（http://www.shanxitv.org/371）。

當前這個例子反映出NTO分析的價值。基于MO討論的話最大的軌道對貢獻只有58.5%，而變換成NTO后最大貢獻達到了86.5%。雖然不完美，不那么接近100%，但還算不錯，起碼有主導的軌道對了，只需討論一對軌道就夠了。不過，對于不少體系的某些激發態，NTO方法并起不到很好效果，即便把MO變換成NTO，單一軌道對的最大貢獻值往往也就70%多甚至更低，光靠這一對軌道還是沒法近似充分反映出電子激發的實際特征。這個時候就必須靠Multiwfn獨家、強大的電子-空穴分析功能了，見《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434），這種分析已經被大量文章采用，如今非常流行。電子-空穴分析比NTO強大得多、普適性強得多，而且使用很方便，因此強烈建議用空穴-電子分析而非NTO！NTO分析唯一的好處僅在于可以給出相位信息，你不需要專門考察相位的話一律應當用空穴-電子分析代替NTO。

這一節的例子是S0到各個單重態激發態之間的NTO分析，對于S0到三重態激發態也可以做NTO分析，只需要在TD里寫上triplet要求計算出來的激發態都是三重態即可。

4 開殼層體系NTO分析實例：苯胺自由基

這里演示一下對一個開殼層體系，苯胺自由基的NTO分析。輸入文件如下
%chk=C:\anilino.chk
# b3lyp/6-31G* TD(nstates=5) IOp(9/40=5)

b3lyp/6-31G* opted

0 2
 C                  0.01820800   -1.80933100    0.00000000
 C                  1.23091700   -1.09944200    0.00000000
 C                  1.22997600    0.28332000    0.00000000
 C                  0.00000000    1.02075200    0.00000000
 C                 -1.22140100    0.26749600    0.00000000
 C                 -1.20379500   -1.11577700    0.00000000
 H                  0.02433600   -2.89560800    0.00000000
 H                  2.17253400   -1.64190800    0.00000000
 H                  2.15204200    0.85613200    0.00000000
 H                 -2.16594900    0.80759600    0.00000000
 H                 -2.13851400   -1.67034400    0.00000000
 N                  0.07314900    2.35943200    0.00000000
 H                 -0.87991900    2.74601300    0.00000000

把chk轉化為fch，載入Multiwfn，進入NTO分析功能，載入Gaussian輸出文件。此例對第5個激發態做分析，輸出信息如下
The sum of square of excitation coefficients:  1.018612
The negative of the sum of square of de-excitation coefficients: -0.018619
The sum of above two values  0.999994
Deviation to expected normalization value (1.0) is  0.000006

The highest 10 eigenvalues of alpha NTO pairs:
   0.540688    0.008875    0.006244    0.000382    0.000273
   0.000245    0.000189    0.000161    0.000100    0.000090
Sum of all alpha eigenvalues:  0.557249

The highest 10 eigenvalues of beta NTO pairs:
   0.390181    0.066067    0.003006    0.000555    0.000330
   0.000291    0.000238    0.000150    0.000127    0.000077
Sum of all beta eigenvalues:  0.461021

Sum of all alpha and beta NTO eigenvalues:  1.018270
 
可見，程序對Alpha和Beta部分分別分析，得到了Alpha NTO和Beta NTO。alpha當中對當前躍遷貢獻最大的一對NTO的貢獻量是0.540688*100%=54.1%，beta當中對當前躍遷貢獻最大的一對NTO的貢獻量是0.390181*100%=39.0%。我們可以讓程序導出fch然后觀看相應的軌道。用Multiwfn載入新產生的.fch文件，進入主功能0，選Orbital info. - Show up to LUMO+10，從文本窗口的信息中很容易就能找到要考察的NTO編號
...[略]
Orb:    24 Ene(au/eV):     0.008875       0.2415 Occ: 1.000000 Type: A
Orb:    25 Ene(au/eV):     0.540688      14.7129 Occ: 1.000000 Type: A
Orb:    26 Ene(au/eV):     0.540688      14.7129 Occ: 0.000000 Type: A
Orb:    27 Ene(au/eV):     0.008875       0.2415 Occ: 0.000000 Type: A
...[略]
Orb:   140 Ene(au/eV):     0.066067       1.7978 Occ: 1.000000 Type: B
Orb:   141 Ene(au/eV):     0.390181      10.6174 Occ: 1.000000 Type: B
Orb:   142 Ene(au/eV):     0.390181      10.6174 Occ: 0.000000 Type: B
Orb:   143 Ene(au/eV):     0.066067       1.7978 Occ: 0.000000 Type: B
...[略]
在軌道觀看窗口里分別選25和26，從而看到的alpha NTO25 -> alpha NTO26的躍遷就是對電子激發貢獻了54.1%的躍遷。而在軌道觀看窗口里分別選141和142，窗口右下角的文本框內容自動就變成了-24和-25，因此這倆軌道對應的是beta NTO24和beta NTO25，前者向后者的躍遷對電子激發貢獻了39.0%。圖示如下：

  
 

到底對此體系用NTO分析比起直接用MO分析有沒有額外優勢？不妨看看MO躍遷的貢獻情況。當前電子激發的MO躍遷情況是如下這樣的（進入主功能18的子功能1 hole-electron分析模塊里的選項10，然后選-2，然后輸入閾值0.05，就會看到組態系數絕對值大于0.05的MO躍遷及其貢獻）
 19A ->    32A   Coeff.:    0.0638   Contri.:    0.4076%
 22A ->    27A   Coeff.:    0.0688   Contri.:    0.4733%
 24A ->    26A   Coeff.:    0.6758   Contri.:   45.6706%
 25A ->    26A   Coeff.:   -0.2878   Contri.:    8.2840%
 25A ->    32A   Coeff.:   -0.0625   Contri.:    0.3904%
 19B ->    25B   Coeff.:   -0.1598   Contri.:    2.5552%
 19B ->    32B   Coeff.:   -0.0712   Contri.:    0.5077%
 22B ->    25B   Coeff.:    0.1833   Contri.:    3.3588%
 22B ->    27B   Coeff.:   -0.0854   Contri.:    0.7293%
 24B ->    25B   Coeff.:   -0.0919   Contri.:    0.8438%
 24B ->    26B   Coeff.:    0.6116   Contri.:   37.4042%
 24A <-    26A   Coeff.:    0.0822   Contri.:   -0.6752%
 24B <-    26B   Coeff.:    0.0721   Contri.:   -0.5197%
可見，如果用MO分析的話，alpha中貢獻第二大的也有8.3%，并沒有小到可以忽略，而轉化成NTO后，alpha當中貢獻第二大的NTO躍遷僅有0.9%了，完全可以忽略。因此NTO對開殼層體系的電子激發分析一樣是有價值的。不過，對開殼層體系需要分別考察Alpha和Beta部分躍遷，總共涉及四個軌道，而用Multiwfn的hole-electron分析的話，就只需要考察hole和electron分布即可，不區分自旋，討論更為省事。

5 做批量NTO分析

經常我們要對一個體系的一大批激發態做NTO分析，雖然也可以比如在Gaussian里用一大堆--link1--分割來做一系列計算得到不同激發態的NTO分析，但是又麻煩又耗時。通過寫shell腳本，可以很方便、快速地讓Multiwfn對一批激發態都做NTO分析、產生記錄NTO軌道的.fch/.molden文件。這里以N-苯基吡咯體系來演示。Multiwfn壓縮包自帶的examples\N-phenylpyrrole.fch和examples\N-phenylpyrrole_ext.out是可以用于NTO分析的輸入文件，此任務用TDDFT計算了5個激發態，我們要一次性對這5個態都做NTO分析，用以下Linux的bash shell腳本可以實現（利用DOS批處理也可以實現相同目的，請自行嘗試）：

#!/bin/bash
cat << EOF > allNTO.txt
18
6
examples/N-phenylpyrrole_ext.out
EOF
for ((i=1;i<=5;i=i+1))
do
cat << EOF >> allNTO.txt
$i
2
S$i.fch
6
EOF
done
./Multiwfn examples/N-phenylpyrrole.fch < allNTO.txt
rm ./allNTO.txt
如果不對此腳本做任何修改，則運行方法是：把此腳本內容復制到文本文件里并命名為allNTO.sh，用chmod 777 allNTO.sh給其加上可執行權限，然后放到Multiwfn的目錄下。進入Multiwfn目錄，運行./allNTO.sh執行此腳本。由于體系不大，一瞬間就會執行完，當你看到當前目錄下產生了S1.fch、S2.fch、S3.fch、S4.fch和S5.fch，說明腳本已執行成功。之后就可以打開相應的fch查看NTO軌道和本征值了。

稍有shell腳本編寫常識的人都很容易理解這腳本是怎么工作的。首先此腳本會產生個allNTO.txt文件，里面包含了在Multiwfn里面要依次敲入的所有命令，然后通過重定向方式調用Multiwfn。i=1;i<=5對應從1到5做循環，讓Multiwfn載入第i號躍遷信息后就產生記錄相應NTO的.fch文件。最后把這個臨時的allNTO.txt刪掉。根據實際情況，自行改一下里面的文件名和循環的上下限就可以用于自己的體系的研究。通過腳本非常靈活方便地使用Multiwfn在《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）里有非常系統的介紹，強烈建議閱讀！

6 關于NTO本征值大于1的問題

時常有人問，他算出來的NTO本征值大于1，結果是否正常、能不能用。這要看具體情況。對于CIS和TDA近似的TDDFT，NTO本征值的范圍一定在[0,1]區間內。而對于TDHF、TDDFT，由于還涉及去激發組態（在NTO方法原文里并沒有明確考慮），NTO本征值有可能大于1。如果只是輕微大于1，比如1.02，就當成1.0就行了。但如果顯著大于1，比如1.5，意味著去激發組態的貢獻巨大，這時候且不說NTO的結果是否合理，可能TDDFT計算本身就是不合理的。這種時候建議對參考態波函數做波函數穩定性測試，確保波函數是穩定的，否則電子激發計算無意義；如果能確保，則建議改用TDA近似的TDDFT。