躍遷密度分析方法-自然躍遷軌道(NTO)簡介

注1：此文發布后Multiwfn也支持了NTO分析，用Multiwfn做NTO比此文的方法方便得多得多得多！，故強烈建議用Multiwfn做NTO分析，而一定不要用Gaussian自己的NTO功能。尤其是處理大批激發態的時候更是必用Multiwfn！用Multiwfn做NTO的方法見《使用Multiwfn做自然躍遷軌道(NTO)分析》（http://www.shanxitv.org/377）。

注2：Multiwfn支持的空穴-電子分析遠比NTO強大得多、普適性強得多，如今非常流行，詳細介紹見《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434）。筆者強烈建議使用空穴-電子分析代替NTO！

躍遷密度分析方法-自然躍遷軌道(NTO)簡介

文/Sobereva @北京科音
First release: 2011-Jun-17   Last update: 2014-Dec-23

1 前言

雖然電子躍遷問題本質是體系電子態與電子態之間的躍遷，但是人們往往喜歡通過軌道模型來描述，也就是把電子躍遷問題以某占據軌道的電子向某虛軌道躍遷的形式來描述，以使躍遷模式更易于考察和理解。使用CIS、TDHF、TDDFT計算激發態時，躍遷模式都可以表達為基態占據軌道與虛軌道間的不同方式躍遷的線性組合，很多情況下，只有一種躍遷所占成分較大，例如Gaussian輸出的一例（RCIS/6-31+G*，丙烯醛）
 Excited State   2:      Singlet-A'     7.0919 eV  174.83 nm  f=0.7713  <S**2>=0.000
      15 -> 16         0.66677
      15 -> 25        -0.13319
這表明，S0->S2這種躍遷方式有88.9%（0.66677^2*2*100%）的程度可以描述為占據軌道15的電子向虛軌道16的躍遷。因此，可以通過分析這種占主導地位的軌道躍遷模式來闡述S0->S2躍遷的本質。

不過，有些情況很難用單一的軌道躍遷模式來描述，例如上面那個體系的S0->S1
 Excited State   1:      Singlet-A"     4.6966 eV  263.99 nm  f=0.0003  <S**2>=0.000
      14 -> 16         0.56407
      14 -> 20        -0.21920
      14 -> 25        -0.24192
      14 -> 29        -0.17712
      14 -> 33        -0.13568
系數最大的是14 -> 16，但它只能描述S0->S1躍遷的63.6%的情況，這不夠充分。然而如果將其它系數不太小的軌道躍遷模式也一起考慮，分析起來會很困難，不便于指認。自然躍遷軌道(Natural transition orbitals, NTO，原文見JCP,118,4775)可以解決這個困難，它通過對分子軌道的變換，使得需要涉及多個軌道躍遷模式才能描述的躍遷問題只用一對兒占據軌道到虛軌道的躍遷就能很好地描述，分析起來容易許多。

2 躍遷密度

先簡要回顧一下躍遷密度的概念。
基態ΨG到激發態ΨE的躍遷密度（無自旋形式，后同）寫為
T(r|r')=T(r1|r1')=N*∫∫...∫|ΨG(x1,x2...xN)><ΨE(x1',x2...xN)| ds1 dx2 dx3... dxN  其中N是總電子數，x是自旋+空間坐標，s是自旋坐標，r是空間坐標。
以基態的分子軌道{ψ}為基的躍遷密度矩陣的i,n矩陣元可寫為
T_i,n=∫∫T(r|r')*ψi`(r)*ψn(r') dr dr'
這里i代表占據軌道標號，n代表虛軌道標號，`是取共軛。T是n_occ*n_vir維矩陣，n_occ和n_vir分別是占據和虛軌道的數目，假設n_vir>n_occ。

PS:這里并不把T寫為n*n維矩陣，n代表全部軌道數。因為這里構成ΨG的各行列式是由ΨE的單行列式通過激發電子得到的，根據Slater-Condon定理，T(r|r')只由各個占據軌道{ψi(r)}和各個虛軌道波函數{ψn`(r')}乘積所構成，T_i,n實際上為ΨE中由i向n激發的行列式的系數乘以N。若令j和m分別代表另外任意一個占據軌道和虛軌道，則T_i,j=0、T_m,n=0，也就是只有n_occ*n_vir的子矩陣不都為0。

對于單粒子算符∑[l]h_l，其期望值<ΨE|∑[l]h_l|ΨG>=∑[i]∑[n]T_i,n*<ψn|h|ψi>。

3 NTO的原理

首先回憶多組態波函數，它是由基態單行列式與對它進行各種激發獲得的其它行列式組合而成的，由于組態數眾多而很難考察。若將密度矩陣對角化，得到的本征向量就是自然軌道，通常只需要本征值最大的n_occ個軌道就能很充分地描述體系密度矩陣，冗余的信息就被去掉了。NTO的思路和做法在某種意義上與此相似，使一大堆軌道躍遷模式轉化為一個“緊湊”的軌道躍遷模式。

首先將TT`矩陣進行酉變換對角化（`仍是取共軛）
U`TT`U=A
U、TT`、A都是n_occ*n_occ矩陣，U是變換矩陣。第i個占據的NTO φi與原先的占據MO {ψj}關系為
φi=∑[j=1→n_occ]U_j,i*ψj
A_i,i就是φi的本征值

然后將T`T進行酉變換對角化
V`T`TV=B
V、T`T、B都是n_vir*n_vir矩陣，V是變換矩陣。第n個虛NTO φn與原先的虛MO {ψm}關系為
φn=∑[m=1→n_occ]V_m,n*ψm
B_n,n就是φn的本征值

這樣，就可以將原先n_occ個占據軌道變換成數目相同的占據NTO，n_vir個虛軌道就被變換成了數目相同的虛NTO。占據以及虛NTO的本征值都小于等于1且大于等于0。每個占據NTO都有一個本征值相同的虛NTO與之對應，因此總共形成n_occ對兒占據→虛NTO躍遷模式。剩下的n_vir減n_occ個虛NTO的本征值都等于0。假設某NTO躍遷模式的本征為0.95，就說明這種NTO躍遷模式可以描述95%的體系電子躍遷模式，或者說這種NTO躍遷模式的組合系數為√0.95=0.975。一般本征值最大的NTO躍遷模式本征值都比較接近于1（除非一些含簡并的高對稱性體系），故一般只需要這一個NTO躍遷模式就足夠展現體系當前電子躍遷的本質特點了。對于CIS，全部占據NTO（或全部虛NTO）的本征值加和為1，代表全部NTO躍遷模式一起可以100%地描述當前的電子躍遷，但對于TDHF/TDDFT由于存在去激發，可能稍微偏離1。

如果用占據的和虛NTO構建躍遷密度矩陣，則可以將躍遷密度矩陣寫成對角的形式，故單粒子躍遷屬性就可以寫為每個NTO躍遷模式貢獻和。其中占主導的NTO躍遷模式的貢獻就可以用于近似解釋躍遷屬性。比如體系的躍遷偶極矩，可以近似用本征值最大的占據NTO到相應虛NTO的躍遷偶極矩來描述。

4 實例

在Gaussian09中已支持NTO分析。這里我們將第一節涉及的丙烯醛S0->S1躍遷做NTO分析。Route section內容如下
# cis/6-31+G* pop(saveNTO,NTO) density=transition=1

density=transition=1代表將基態到第一激發態的躍遷密度矩陣傳遞給L601模塊（波函數分析模塊）用于分析，pop里NTO代表對這個傳來的S0->S1躍遷密度矩陣生成NTO，saveNTO代表將生成的所有NTO都保存到checkpoint文件里，就可以用gview，或者轉化為fch文件后用Multiwfn (http://www.shanxitv.org/multiwfn) 觀看NTO軌道了，而且還可以用Multiwfn做軌道成分分析。

Gaussian里對NTO的排序方式是：對于占據NTO，本征值從低往高排，對于虛NTO，本征值從高往低排。下面是輸出信息的一部分，占據軌道部分最后一個（第15號）是0.99276，這就是本征值最大的占據NTO，虛軌道部分第一個（第16號）就是本征值最大的虛NTO。可見每個占據NTO都有一個本征值相同的虛NTO相對應。0.99276這個值幾乎達到了上限，它說明S0->S1躍遷的高達99.2%的內涵都可以只用15號軌道到16號NTO的躍遷來描述，其它的NTO躍遷模式，如14到17號（占0.6%）、13到18號（占0.08%）等等都可以忽略。

 Alpha  occ. eigenvalues --    0.00000   0.00000   0.00000   0.00000   0.00000
 Alpha  occ. eigenvalues --    0.00000   0.00000   0.00001   0.00003   0.00005
 Alpha  occ. eigenvalues --    0.00008   0.00018   0.00077   0.00611   0.99276
 Alpha virt. eigenvalues --    0.99276   0.00611   0.00077   0.00018   0.00008
 Alpha virt. eigenvalues --    0.00005   0.00003   0.00001   0.00000   0.00000
 Alpha virt. eigenvalues --    0.00000   0.00000   0.00000   0.00000   0.00000
...略

把chk轉換為fch文件，再用Multiwfn載入之，進入主功能0查看分子軌道，第15和第16號NTO如下所示

我們通過Multiwfn基于Becke方法來做一下15號NTO的軌道成分分析。關閉觀看分子軌道的界面，然后輸入
8   //軌道成分分析
9   //Becke方法
15   //15號軌道
我們看到第15號NTO中氧的成分達到了70%，因此可以指認S0->S1躍遷是氧的孤對電子向π反鍵軌道躍遷。用NTO的時候考察一對兒軌道的躍遷就說明問題了，而原先使用分子軌道描述時需要同時考慮好幾對兒軌道的躍遷，明顯NTO使躍遷分析方便多了。

如果接下來想對比如S0->S2做NTO分析，當然可以用cis/6-31+G* pop(saveNTO,NTO) density=transition=2，但是這要完全地重算一遍，比較費時。較好的辦法是寫cis(read)/6-31+G* guess=read pop(saveNTO,NTO) density=transition=2，這說明在SCF過程中直接讀取check文件里已經收斂的SCF波函數，而在CIS過程中（需要Davidson迭代求解）也直接讀取check文件里的初猜，這樣做會比起重算一遍速度會快不少。

5 其它

NTO對于一些情況，特別是簡并的情況，效果并不好，也就是多個NTO對兒的本征值都不小，沒有哪對兒NTO是主導的從而可以較好描述電子的躍遷，此時仍需要通過兩個或更多NTO對兒才能描述電子躍遷，相當于NTO分析失去意義了。在Multiwfn程序(http://www.shanxitv.org/multiwfn)的主功能18里的子功能1可以做電子-空穴(electron-hole)分析，會給出電子以及空穴的分布圖形。假設躍遷可以被一對兒NTO所完美描述，那么這種電子-空穴分析與NTO分析給出的結論是一致的，電子分布和空穴分布分別相當于那兩個NTO軌道波函數的平方。然而當單個NTO對兒不足以定性描述躍遷時，就必須用電子-空穴分析了，也就是說Multiwfn的電子-空穴分析是一種比NTO更普適的分析躍遷的方法，這種分析方法詳見手冊3.21.1節的理論介紹和4.18.1節的示例。另外，如果對于一個體系你需要分析多個激發的話，Multiwfn的這個功能比起NTO用起來也更方便，因為只需要用Gaussian計算一次即可，然后在Multiwfn里可以直接選擇要分析哪個激發態，而不用像NTO一樣每次還得重新用Gaussian算一遍。

值得一提的是Mayer在CPL,437,284里提出的方法和NTO在出發點和目的上如出一轍，只是推導過程不同而已，但實際上也不難證明等價性。Mayer的推導處理的不是躍遷密度，而是CI系數。對應的程序CIS-T可以從這里下載http://occam.chemres.hu/programs/，不過也只能處理Gaussian的輸出，由于G09已經內建了NTO，CIS-T就用處不大了。

NTO同樣可以分析開殼層體系的躍遷，意義與上面介紹的是相同的，用前文處理無自旋的躍遷密度矩陣的方法分別處理alpha-alpha躍遷密度矩陣和beta-beta躍遷密度矩陣即可。

補充說明：很多人在NTO和分子軌道的序號方面沒搞明白。這里再強調一下：NTO和MO的序號沒有任何直接對應關系！如文中所說，Gaussian里NTO是根據相應的TT`、T`T矩陣的本征值排序的，而MO是按照能量從低到高排序的，序號規則完全不同。而且NTO是從MO變換而來（或者說每個NTO是由諸多MO混合得到），形狀上也沒法直接對應。一種電子躍遷模式多數情況都可以僅靠本征值最高的占據NTO向相應本征值的虛NTO躍遷來近似描述。例如對于本文討論的丙烯醛（15個占據軌道），不管是哪種躍遷，無論是S0->S1也好還是S0->S2也好或是向更高的激發態躍遷也好，按照Gaussian對NTO的排序方式，總是能近似視為是NTO 15->NTO 16的躍遷。然而，這些躍遷用MO的躍遷來表示就很不相同了，比如S0->S2主要以MO 15->MO 16來表示（HOMO->LUMO），而占S0->S1最大成分的則是MO 14->MO 16。另外，NTO也沒有軌道能量的概念，絕對不能由于電子激發模式可以被NTO 15->NTO 16的躍遷所很好地表達就說成是HOMO->LUMO的躍遷，這是大錯特錯。