使用Multiwfn繪制躍遷密度矩陣和電荷轉移矩陣考察電子激發特征（含視頻）
本文的實例操作部分有視頻演示，見：https://www.bilibili.com/video/av35812201/
使用Multiwfn繪制躍遷密度矩陣和電荷轉移矩陣考察電子激發特征
文/Sobereva @北京科音
First release: 2018-Nov-12 Last update: 2021-May-27
0 前言
躍遷矩陣（transition matrix）是一個泛稱，表示各種蘊含了電子躍遷特征信息的矩陣。其中最常被研究的是躍遷密度矩陣（transition density matrix, TDM），其它的各種躍遷矩陣都是由TDM衍生而來，或者與TDM存在密切關系。大多數躍遷矩陣有兩種表現形式，一種是實空間形式，可以通過等值面圖、平面圖來表現；另一種是一般意義的矩陣形式，矩陣元的標號可以對應于基函數序號、原子序號、片段序號，可以通過繪制熱圖（即填色矩陣圖）來直觀地展現。通過這些圖可以一目了然地了解體系各個位點和位點間的耦合對電子激發的影響，對于幫助了解電子激發內在特征十分有幫助。本文介紹的研究方法和Multiwfn中的空穴-電子分析在原理上存在一定共通性，在展現形式上也有互補性，可以相結合討論來更好地表征電子激發特征。
本文的主要目的是介紹怎么用Multiwfn產生和分析各種類型躍遷矩陣的等值面圖和熱圖，使得讀者可以在研究電子激發問題時擁有更強大的手段。Multiwfn可在其主頁http://www.shanxitv.org/multiwfn免費下載。如果對Multiwfn不了解，建議參看《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）。強烈建議在閱讀本文前先仔細閱讀《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434）和《在Multiwfn中通過IFCT方法計算電子激發過程中任意片段間的電子轉移量》（http://www.shanxitv.org/433），這十分有助于更好地理解本文的內容。如果讀者對于電子激發計算比較陌生，務必先看《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）。
值得一提的是很久以前筆者專門寫過一篇相應博文《繪制躍遷密度矩陣平面圖分析電子躍遷》（http://www.shanxitv.org/136），這篇文章介紹的是如何在Multiwfn中繪制原子TDM的熱圖，此文發布后有大量研究文章都使用Multiwfn做了這種分析，比如J. Phys. Chem. C, 121, 2574 (2017)、RSC Adv., 7, 19576 (2017)、J. Mol. Model., 23, 28 (2017)等等。但現在來看這篇文章相對于目前版本的Multiwfn來說已經完全過時了，而且寫得很糙，有些敘述也不太準確，所以看本文后就不要再看那篇文章了。
在講操作和實例之前，先介紹一下躍遷密度矩陣(TDM)及相關概念。
1 基礎知識
密度矩陣是對特定的某個電子態而言的，它是這個電子態波函數的一種變相表現形式。我們一般說的密度矩陣特指一階約化密度矩陣，通過它可以獲得這個態的各種單電子性質，諸如電子分布、偶極矩、靜電勢等。TDM則是對兩個態之間的躍遷而言的，基于它可以獲得躍遷過程中的一切單電子性質，諸如原子躍遷電荷、躍遷電/磁偶極矩等。TDM有兩種表現形式，一種是實空間形式，另一種是一般意義上的矩陣形式（基函數表象下的形式），這兩種表現方式在討論實際問題時各有所長，下面依次介紹。類似地，TDM所衍生的躍遷偶極矩密度矩陣也有這兩種表現形式，它與躍遷偶極矩密切相關，下面也一起介紹。
注：不考慮旋軌耦合時，TDM既有空間部分也有自旋部分，不同自旋多重度的電子態之間的TDM精確為0，但本文說的TDM及相關的矩陣，以及在Multiwfn中所產生的，一律是空間部分，因此即便是比如單重態躍遷到三重態對應的TDM也不是零矩陣。
1.1 TDM的實空間形式
對多電子體系，基態與激發態間的TDM的實空間形式T(r;r')通過下圖第一個式子得到。
式中Φ0是基態波函數，Ψexc是某個激發態波函數，x是電子的自旋+空間坐標，σ是自旋坐標，r是空間坐標。之所以叫做“實空間形式”，是因為這樣的TDM的兩個指標是三維空間的坐標矢量。
如果計算激發態用的方法是單參考方法，比如常見的TDDFT、CIS、ZINDO，則激發態波函數通過各種單激發組態函數的線性組合來描述，此時T(r;r')可以更確切地寫為上圖第二個式子，其中a和i分別是空軌道和占據軌道序號，w是組態函數的系數。
由于r含有三個分量，因此T(r;r')是六維函數，沒法簡單地通過圖像表示。如果我們取TDM的對角元，即讓r=r'，此時T(r;r')就變成了T(r)，即上圖的第三個式子，這叫做“躍遷密度”。由于它是三維函數，可以很容易地作圖來考察。
躍遷密度看起來有點抽象，為了便于理解，我們可以考慮一個最簡單的情況，即電子激發恰可以完美地通過一個占據軌道i向一個空軌道a的躍遷來描述，此時T(r)=φi(r)*φa(r)。顯然躍遷密度有的地方為正有的地方為負，但正負不是我們這里關心的，這里我們只關心大小。躍遷密度越大的地方，對應i、a軌道在此處重疊越顯著，即這兩個軌道波函數的乘積的絕對值大；而i、a重疊小的地方，顯然躍遷密度也會很小。實際的電子激發一般都不可能被一對軌道躍遷所完美描述，而應當廣義地描述為“空穴→電子”。推廣一下，我們也可以近似認為躍遷密度大的地方正是空穴和電子重疊程度大的地方。躍遷密度小的地方并不是說這樣的地方在電子激發過程中沒有牽扯到，只不過是說空穴和電子并沒有“同時”在這個地方有較大分布。假設一個電子激發是A區域向B區域的電子轉移激發，那么看到的躍遷密度圖應該是在A和B相交的地方有較大分布，而在其余地方則較小。值得一提的是，由于分子軌道間的正交性，T(r)在全空間積分必定恰為0。
將躍遷密度分別與x,y,z坐標變量相乘，就得到了躍遷偶極矩密度的x,y,z三個分量，如下所示。
比如對Tx(r)進行全空間積分，就恰等于躍遷偶極矩的X分量。因此，利用躍遷偶極矩密度，可以十分方便地考察體系各個位置對躍遷偶極矩的貢獻。眾所周知，決定了電子態間躍遷概率的振子強度正比于躍遷偶極矩的模的平方（這點在http://www.shanxitv.org/314中介紹了），因此對躍遷偶極矩密度作圖，就可以一目了然地了解體系哪些區域對兩個態間的躍遷概率有什么樣的貢獻。
上面所說的躍遷偶極矩特指躍遷電偶極矩，這也是我們平時最關心的躍遷偶極矩的形式。實際上還有躍遷速度偶極矩、躍遷磁偶極矩等。Multiwfn也可以給出躍遷磁偶極矩密度，公式見Multiwfn手冊3.21.1.2節，將之全空間積分就是躍遷磁偶極矩了，這和用于繪制ECD光譜的轉子強度直接相關。
1.2 TDM的常規矩陣形式
分子軌道本身是三維實空間函數，但是定義了基函數后，就可以將之通過各個基函數的展開系數來表示。類似地，在一套特定基函數下，TDM的兩個指標就不再是連續的實空間坐標r和r'了，而是基函數的序號。此時基態與激發態之間的TDM可以通過下圖第一個式子計算，它與實空間形式的TDM間可以通過下圖第二個式子轉化。
通過sum-over-states (SOS)方法計算體系的一些屬性等場合需要利用到激發態之間的躍遷信息，因此激發態之間的TDM也是有實際價值的，Multiwfn也能給出，公式這里就不多說了。
一種常用的直觀展現各個矩陣元大小的方式是熱圖(heat map)，即圖像里每個格子對應一個矩陣元，根據矩陣元的數值對格子進行著色，這樣一目了然就知道哪些矩陣元數值較大，以及是正是負，TDM也可以以這種方式展現。
TDM看起來比較抽象，如何通過觀看它的熱圖來提取有意義的信息？我們可以看一個最簡單的情況，假設體系只有兩個基函數，而且電子激發可以完美地用軌道i→軌道a的躍遷來描述，此時的TDM如下，其中比如C_2a就代表a軌道中第2個基函數的展開系數。注意為了和文獻中繪制的TDM的熱圖的習俗相一致，此處矩陣元的序號順序和一般數學上的矩陣序號順序不同。
我們先看TDM的對角元。根據上面的表達式可見，如果(1,1)矩陣元的絕對值很大，必然1號基函數同時在i和a中的系數很大，可以廣義地說1號基函數對空穴和電子都有很大貢獻，或者說空穴和電子在1號基函數上有顯著重疊。雖然密度矩陣是對稱矩陣，但TDM一般不是對稱矩陣。根據上式可見，如果(1,2)矩陣元比較大，則說明1號基函數和2號基函數分別在i和a中的展開系數的絕對值較大，可以廣義地說1號基函數參與空穴較大而且與此同時2號基函數參與電子較大，根據筆者提出的IFCT分析的思想，從物理上可以理解為這種電子激發會導致體系的電子從1號基函數向2號基函數轉移。而如果(2,1)矩陣元比較大，則說明電子激發會造成體系的電子從2號基函數向1號基函數轉移。
相對于上一節提到的躍遷密度圖，顯然TDM熱圖還蘊含了更多信息，即它的非對角元還體現了電子轉移特征，這在躍遷密度圖上不直接體現。
1.3 原子TDM、片段TDM
上述這種TDM對應的熱圖并不是很好分析，因為基函數只有數學意義，通常不是我們感興趣的層面。我們感興趣的通常是電子激發牽扯到體系中的哪些原子，以及電子是怎么在原子間轉移的。為此，我們可以把TDM人為地“收縮”成原子TDM，之后矩陣元的序號就對應原子序號了。收縮方式并不唯一，Multiwfn支持下面幾種方法，從形式上看大同小異，被不同文獻所使用：
其中μ、ν是基函數序號，A、B是原子序號。TDM矩陣元有正有負，但正負號沒太大物理意義，還會導致直接加和時相互抵消，因此上述把TDM收縮成原子的方法都是先通過取絕對值或者求平方把TDM矩陣元的符號去除后再加和成原子。由于氫原子通常對化學上感興趣的電子激發貢獻極小，因此為了矩陣的熱圖看著緊湊，通常在作圖時忽略氫原子。
注：按照上述方法2,3,4產生的p矩陣是真正意義上的原子TDM，而第一種方法產生的p矩陣，其實嚴格來說不叫原子TDM，它對應的是所謂的correlated electron-hole probability diagram (CEHPD)這種圖的矩陣，(A,B)矩陣元被解釋為同時在A和B原子上發現空穴和電子的概率，見比如J. Chem. Phys., 113, 10002 (2000)和J. Am. Chem. Soc., 129, 14257 (2007)的例子。但不管上述哪種方法構建的p矩陣，在定性特征上都是相似的，繪制出的熱圖的效果一般也定性相同。為了敘述方便，下面都用“原子TDM”來統一稱呼。
為了討論方便，還可以構建并繪制片段TDM，它通過把原子TDM矩陣元直接加和來構建。
原子/片段TDM的一般形式可以這么示意：
類似于對TDM的分析，我們對此矩陣可以近似這么理解：如果對角元(I,I)大，就說明電子和空穴同時在I位點上具有較大分布；如果非對角元(I,J)大，就說明空穴在I位點上大，同時電子在J位點上大，因此可以認為此電子激發伴隨著I向J的電子轉移。如果矩陣的左上部分和右下部分是基本對稱的，即(I,J)≈(J,I)，就暗示電子激發時I→J和J→I的轉移量基本相同，基本沒有出現電荷的凈轉移。
值得一提的是，一些文獻使用的TDM是按照下式對稱化后的形式
我不建議效仿這些文獻用對稱化后的TDM討論，因為基于這樣的TDM構建的原子/片段的躍遷矩陣的對角元就不再體現電子轉移方向的信息。此時的(A,B)或(B,A)矩陣元大，只能簡單地說A-B之間存在電子轉移，這在一些文章中也被稱為A和B存在相干(coherence)。
1.4 原子/片段的電荷轉移矩陣
要注意，上面討論TDM時提到的“電子”和“空穴”只是一種抽象的概念，這和Multiwfn中的空穴-電子分析方法中所定義的空穴和電子雖然在物理意義上類似，但絕對不是一碼事，因為Multiwfn的空穴-電子分析里的空穴和電子是通過數學形式嚴格定義的。也因此，TDM圖所展現的特征，和Multiwfn產生的空穴、電子的分布特征，雖然大多數情況下能對應起來，但并非總是能很好對應。
在筆者來看，比原子/片段TDM物理意義明顯更明確，也同時更有實際價值的是原子-原子或片段-片段的電荷轉移矩陣(charge transfer matrix)。電荷轉移矩陣雖然和TDM的關系不是特別近，但也可以算是原子躍遷矩陣的一種。這種矩陣是筆者提出的IFCT分析方法的副產物，它的定義很簡單，其矩陣元(I,J)=Θ(I,hole)*Θ(J,electron)。這里比如Θ(I,hole)就是指I原子或片段在hole中所占比率。比率的計算方法并不唯一，Multiwfn也支持不止一種算法，詳見前述的《在Multiwfn中通過IFCT方法計算電子激發過程中任意片段間的電子轉移量》一文。
根據IFCT方法的思想，(I,I)矩陣元直接對應電子激發時在I位點發生的電子重分布量，而(I,J)矩陣元直接對應從I位點向J位點轉移的電子量。電荷轉移矩陣顯然也可以通過熱圖表現，這種分析方法在目前撰文時筆者還沒發表，但筆者認為，有了電荷轉移矩陣的熱圖，其實就沒太大必要再用原子/片段TDM熱圖分析了，因為電荷轉移矩陣比原子/片段TDM物理意義明顯更明確，解釋起來更清楚，而且和Multiwfn給出的空穴、電子分布圖形總是精確對應。
1.5 原子躍遷偶極矩矩陣
躍遷偶極矩矩陣的矩陣元定義為對應的TDM矩陣元與相應基函數間的偶極矩積分的乘積。躍遷偶極矩矩陣有X,Y,Z三個分量，諸如X分量矩陣的所有矩陣元加恰為體系躍遷偶極矩的X分量。也可以將躍遷偶極矩矩陣收縮成基于原子的形式，原子躍遷偶極矩矩陣的X,Y,Z分量定義為
諸如X分量矩陣的(A,A)矩陣元就是完全由A原子自己對體系躍遷偶極矩X分量的貢獻，而(A,B)則對應于A、B原子產生的聯合貢獻。還可以定義原子總躍遷偶極矩矩陣，其中每個矩陣元是相應X,Y,Z分量的平方和。原子躍遷偶極矩矩陣還可以進一步收縮成基于片段的形式。通過繪制原子/片段躍遷偶極矩矩陣的熱圖，對于解釋躍遷偶極矩的大小、分析主要來源是有益的。
2 Multiwfn中與躍遷密度相關的分析
2.1 相關功能
在Multiwfn中，與本文主題有關的功能如下，都在主功能18里：
(a)子功能1：這是空穴-電子分析模塊，在計算空穴和電子格點數據的同時會順帶計算出前面提到的躍遷密度和躍遷偶極矩密度格點數據，可以直接觀看其等值面、導出cube文件。
(b)子功能8：這是IFCT分析模塊，可以導出前述的原子-原子電荷轉移矩陣到當前目錄下的atmCTmat.txt。
(c)子功能9：產生基態與激發態，以及激發態與激發態之間的TDM，可以導出為當前目錄下的tdmat.txt。也可以導出TDM.fch，此文件中記錄密度矩陣的段落對應于新產生的TDM。
(d)子功能11：可以計算各個原子對基態與激發態間的躍遷偶極矩的貢獻（《使用Multiwfn+VMD繪制片段貢獻的躍遷偶極矩矢量》http://www.shanxitv.org/396一文用到了此功能），還可以導出原子躍遷偶極矩矩陣到當前目錄下的以AAtrdip開頭的.txt文件中。
(e)子功能2：這是本文的重點，用來繪制原子/片段躍遷矩陣的熱圖。躍遷矩陣可以以不同方式得到：可以讓此功能直接生成基態到激發態的TDM，也可以從Gaussian輸出文件中讀取基態到激發態的TDM，也可以從前述的tdmat.txt或AAtrdip.txt或atmCTmat.txt中讀取相應類型的矩陣。
以上功能具體使用細節在手冊3.21節的相應章節里有非常詳細易懂的說明，在此就不細說了，后文將通過實例進行演示。
2.2 輸入文件
使用上述主功能18的子功能1、8、9、11，以及使用子功能2時打算直接產生TDM的話，需要以下兩類文件，更具體說明參見Multiwfn手冊3.21節開頭。
·第一類文件：記錄了參考態波函數的文件，是剛啟動Multiwfn時就需要載入的
·第二類文件：記錄了激發態組態系數的文件，是進入相應分析功能時需要載入的
對于Gaussian用戶，一般就用TDDFT算激發態即可（用CIS、TDHF、TDA-DFT亦可）。計算時候應當帶IOp(9/40=4)關鍵詞使得絕對值大于1E-4的組態系數都輸出出來以使得分析結果可靠。把算完得到的chk轉換為fch文件后就可以作為第一類文件，而Gaussian輸出文件可作為第二類文件使用。這些分析絕不僅限于Gaussian用戶能用，諸如GAMESS-US、Firefly、ORCA等用戶也可以用，詳見手冊3.21節開頭的說明。
使用子功能2繪制熱圖時如果從tdmat.txt或AAtrdip.txt或atmCTmat.txt中直接讀取相應類型的矩陣，那么第一類文件同前，載入第二類文件的時候輸入這些txt文件的路徑即可（文件名不能改為別的，否則程序不知道讀的是哪類矩陣）。
使用子功能2繪制熱圖時如果想從Gaussian輸出文件中直接讀取TDM，那么第一類文件同前，載入第二類文件的時候輸入Gaussian輸出文件路徑即可。此時Gaussian計算時必須寫density=transition=x IOp(6/8=3)關鍵詞，代表把基態到第x激發態的TDM輸出到輸出文件里，而IOp(9/40)則不再需要寫。這種用法一般沒必要用，也只會給出對稱化后的TDM。僅當你的體系很大，不得不使用ZINDO半經驗方法計算的時候，才需要以這種方式把ZINDO級別的TDM提供給Multiwfn繪制成熱圖（Multiwfn自己沒法產生ZINDO這些半經驗方法對應的TDM）。
3 實例
本文用的Multiwfn 3.6(dev)是2018-Nov-12更新的版本，絕對不要用更老版本，Gaussian用的是G16 A.03版。
3.1 繪制躍遷密度和躍遷偶極矩密度圖
首先以一個簡單分子N-phenylpyrrole（N-苯基吡咯）為例演示一下怎么繪制和討論躍遷密度和躍遷偶極矩密度圖。下面用到的.fch和.out文件是使用cam-b3lyp/6-31+G(d) TD(nstates=5) IOp(9/40=4)關鍵詞產生的。
啟動Multiwfn然后輸入
examples\excit\N-phenylpyrrole.fch //文件在Multiwfn自帶的examples文件包里
18 //電子激發分析
1 //空穴-電子分析
examples\excit\N-phenylpyrrole.out
1 //考察基態(S0)到第1個激發態(S1)的躍遷
1 //計算和可視化空穴、電子、躍遷密度等函數
2 //中等質量格點
算完之后，屏幕上輸出了很多信息，大部分是與空穴-電子分析相關的，其中有一條和當前主題有一定關系：
Transition dipole moment in X/Y/Z: -0.000021 -0.000045 1.767332 a.u.
這是基于躍遷偶極矩密度格點數據在全空間做的積分值，即基態到第1激發態的躍遷偶極矩，其數值和Gaussian輸出文件里的十分接近，以下是N-phenylpyrrole.out的第773行：
state X Y Z Dip. S. Osc.
1 0.0000 0.0000 1.7813 3.1729 0.3935
二者十分接近體現出當前研究用的格點設定是足夠合理的，因此之后看到的躍遷密度和躍遷偶極矩密度的等值面圖是可以充分說明實際問題的。
在后處理菜單中選選項5顯示出來躍遷密度的等值面圖，如下圖左半邊所示。如果選選項6，再選擇一個分量，比如這里選Z，就會顯示出躍遷偶極矩密度Z分量的等值面圖，如下圖右半邊所示。等值面數值用的是默認的0.001。
圖中綠色和藍色分別對應正值和負值等值面。躍遷密度等值面遍布整個體系，暗示當前研究的電子激發是整體躍遷，空穴和電子都分布在整個體系上，這點也可以通過繪制空穴和電子等值面圖進一步確認（后處理菜單中選3）。躍遷偶極矩密度Z分量圖實際上就是躍遷偶極矩密度乘上Z坐標再取負值后對應的圖（圖中Z軸由下指向上），從圖中可見體系各個地方對躍遷偶極矩Z分量都有貢獻，由于圖中綠色區域顯著大于藍色區域，這解釋了為什么當前的電子激發的躍遷偶極矩Z分量是明顯的正值（1.767 a.u.）。
下圖顯示的是躍遷偶極矩密度X分量的等值面圖，X是垂直于分子平面的方向。可以看到分子平面左右兩側的等值面顏色恰好相反，正負部分精確抵消，這是為什么當前體系的躍遷偶極矩的X分量精確為0。
我們再來看S0到S5的躍遷情況。按照屏幕上的提示用相應選項退出空穴-電子分析模塊，再次進入，讓你輸入要考察的激發態的時候選擇5，然后按照前述步驟繪制躍遷密度、躍遷偶極矩密度以及空穴&電子圖，分別如下所示
從空穴&電子圖上看，此電子激發是明顯的電荷轉移激發，因為空穴幾乎都分布在吡咯上，電子幾乎都分布在苯環上，所以電子和空穴的重疊程度低，這導致體系各處的躍遷密度的數值都不太大，也因此上圖的S0→S5的躍遷密度等值面圖看起來等值面范圍比起S0→S1的要小得多。當前的電子激發的躍遷偶極矩Z分量為1.06 a.u.，為什么是明顯的正值也是從Tz(r)圖上一目了然地就可以了解，主要是因為分子兩端有顯著的正貢獻，明顯超過了其它一些零星區域的負貢獻。
上面的例子展現了躍遷偶極矩密度對于討論體系局部對躍遷偶極矩貢獻的重要價值。但要注意的是，躍遷偶極矩密度是依賴于原點的選擇的，比如體系在Z正方向整體平移10埃，那么得到的躍遷偶極矩密度Z分量圖就會和之前顯著不同，原因從其表達式上很容易理解，因為它是在T(r)的基礎上乘了個-z項（也因此，越接近原點的區域，躍遷偶極矩密度肯定傾向于越小）。而躍遷密度雖然與躍遷偶極矩之間的關系相對間接一些，但好處是其圖像完全不依賴于原點的選取。
在筆者來看，從表象到深層原因，有這么個順序：
吸收/發射強度←振子強度←躍遷偶極矩←躍遷偶極矩密度←躍遷密度←空穴&電子分布
就算我們不在文章里用躍遷(偶極矩)密度圖來討論，但腦子里有躍遷(偶極矩)密度的概念，對于我們想明白一些問題的本質也是極有幫助的。
如果在空穴-電子分析模塊里先選一次選項"-1 Toggle calculating transit. magnetic dip. density in option 1"將之狀態切換為Yes然后再計算格點數據，則算完格點數據后屏幕上會輸出躍遷磁偶極矩矢量，并且可以在后處理菜單中選相應選項繪制躍遷磁偶極矩密度圖。
躍遷密度、躍遷電/磁偶極矩密度也都可以通過后處理菜單中的選項導出cube文件，通過它還可以在Multiwfn里將之繪制成曲線圖、平面圖。過程參考《使用Multiwfn計算（超）極化率密度》（http://www.shanxitv.org/305）中的相應做法。
對大多數實際研究的體系，躍遷偶極矩并不是像此例這樣恰好沖著某個笛卡爾軸的，若我們想通過一張躍遷偶極矩密度圖就能討論它的貢獻來源，應當旋轉體系讓體系躍遷偶極矩平行于某個笛卡爾軸，比如X軸，然后再重新用諸如Gaussian做一次電子激發計算并結合nosymm關鍵詞，之后再用Multiwfn繪制并分析躍遷偶極矩密度X分量圖即可。具體做法見《讓體系(躍遷)偶極矩平行于某個笛卡爾軸的方法》（http://www.shanxitv.org/507）。
3.2 繪制基態到激發態的原子TDM的熱圖
本例我們演示繪制下圖所示的Donor-pi-Acceptor型長鏈分子NH2-C8-NO2的基態到激發態的原子TDM的熱圖，gjf,、out、fchk文件在examples\excit\NH2_C8_NO2目錄下都提供了，Gaussian計算用的關鍵詞為CAM-B3LYP/6-31G* IOp(9/40=4) TD(nstates=10)。一般來說原子TDM的熱圖只適合討論一維鏈狀體系，這樣熱圖坐標軸上的序號才和原子位置容易對應起來，如果不是一維體系的話建議改用下一節示例的片段TDM圖討論。
前面提到，原子TDM的熱圖里一般是忽略掉氫原子的，因此為了讓熱圖里的序號和結構圖里的序號直接對應，在計算前已經把氫原子的序號排到了非氫原子的后頭，做法是在gview里，進入Edit - Atom List，選Edit - Reorder - All atoms: Hydrogens Last，然后再保存輸入文件（如果你計算前沒這么做也沒關系，計算之后再這么在gview里操作一下，非氫原子序號也能和Multiwfn繪制的非氫原子構成的原子TDM熱圖對應上）。
啟動Multiwfn，然后輸入
examples\excit\NH2_C8_NO2\NH2_C8_NO2.fchk
18 //電子激發分析
2 //繪制原子/片段躍遷矩陣的熱圖
examples\excit\NH2_C8_NO2\NH2_C8_NO2.out
1 //考察S0→S1激發。然后程序會產生這個躍遷對應的TDM
n //不對剛產生的TDM做對稱化
1 //用第1種模式（見本文1.3節的說明）將TDM收縮為原子TDM，此做法比較常用效果也好
1 //繪制熱圖
此時看到下圖。默認的色彩刻度下限為0，上限為矩陣元的最大值。
如前所述，這種圖的對角元的大小體現在什么位置上電子和空穴同時有顯著分布，而考察非對角元的時候，要先看橫坐標（對應空穴位置），再看縱坐標（對應電子位置），根據熱圖中相應位置的數值大小，就可以知道電子在體系各位點間是怎么轉移的。此例圖中的對角線大部分都被綠色或紅色包圍，這立刻體現出這個躍遷必然是全局激發，激發的電子會牽扯幾乎整個體系。此圖不是沿對角線左上和右下對稱的，對角線左上部分的數值整體比右下角部分更大，而且主要是在對角線附近比較大，因此此圖反映出電子激發過程中各非氫原子上的電子會往與之近鄰的原子轉移，而且序號小的原子往序號大的原子轉移的程度比序號大的原子往序號小的原子轉移的更大。由于此體系非氫原子的序號順序是從氨基向硝基排過去的，因此可以推測這種S0→S1激發導致電子從氨基端往硝基端整體移動。
如果覺得TDM的熱圖還有些抽象，可以結合空穴-電子分析給出的空穴&電子等值面圖考察，如下所示
確實，從空穴&電子等值面圖看到的情況和基于TDM熱圖做出的結論是基本一致的（盡管從原理上并不直接對應，因為此圖展現的空穴和電子是嚴格定義、嚴格計算的，而分析TDM圖時提及的空穴和電子只是抽象、含糊的概念）。
如果想繪制其它激發態的TDM熱圖就先退出熱圖繪制功能，再次進入時選擇相應的態即可。繪制熱圖的功能的后處理菜單中有很多選項，含義不言自明，手冊里也講了，這里對其中幾個專門提一下。"4 Toggle if taking hydrogens into account"用來切換是否考慮氫原子，如果選一次切換為Yes，再次作圖，看到的圖就是下面這樣。氫原子序號是13~22，確實從圖上看氫原子沒怎么參與電子激發，對應的矩陣元數值都很小，因此把氫納入S0→S1激發的熱圖確實沒有什么意義。
再看熱圖的格子之間插值對結果的影響。默認情況下插值10次，這樣得到的圖像比較平滑。如果選擇6 Set the number of interpolation steps between grids然后輸入1，就相當于不做插值，得到的圖像如下所示。此時每個格子正好對應一個矩陣元，雖然對應關系精確，但看著就沒那么舒服了。默認用的10次插值已經夠大了，設得更大也可以，繪圖時間會變得更長，而比默認設定所得圖像并沒什么明顯改進。
后處理菜單有個選項"Switch if normalizing the sum of all elements to unity"，如果將之狀態切換為"Yes"，則程序會給矩陣數據乘上一個系數，使得所有矩陣元加和數值為1。這個設定對于當前這個電子激發的圖沒什么明顯影響，但對于某些電子激發，不開這個選項的話可能矩陣元最大值也是非常小的，而開了這個選項后可以顯著增大所有矩陣元的數量級，使得熱圖效果較好時色彩刻度上限不需要設得特別小。
我們再來看看構建原子TDM的時候選擇模式2、3、4的時候圖像效果分別各是什么樣，如下所示
從圖像上看，以上三張圖和之前基于模式1構建的原子的TDM的熱圖傳達的信息很類似。模式2、3的圖像效果也都還可以，都可以用，但模式4對應的圖不是很理想，由于其構建對角元和非對角元用的公式不同，導致圖像上過度強調對角元了，而非對角元的特征展現得不是非常清楚。
我們再來繪制一下原子-原子電荷轉移矩陣對應的熱圖，這種圖和空穴&電子等值面圖在原理上是嚴格對應的，可以嚴格對比。啟動Multiwfn，依次輸入
examples\excit\NH2_C8_NO2\NH2_C8_NO2.fchk
18 //電子激發分析
8 //IFCT分析
1 //類Mulliken方式計算空穴、電子的成份
examples\excit\NH2_C8_NO2\NH2_C8_NO2.out
1 //考察S0→S1激發
-1 //導出原子-原子電荷轉移矩陣到當前目錄下的atmCTmat.txt中
2 //繪制原子/片段躍遷矩陣的熱圖
atmCTmat.txt
1 //將剛載入的atmCTmat.txt里的矩陣繪制為熱圖
此時看到的圖像如下
此圖和之前基于模式2給出的原子TDM的熱圖的基本特征較為相似，但也有不可忽視的差異。當前這個圖的每個非對角元比較嚴格地展現了原子間的電子轉移量，我認為比原子TDM圖更有意義。一列一列地看此圖的話，可以直觀地看到碳鏈上的每個原子都往它前端和后端的原子上轉移了電子，往硝基一側轉移的明顯更多一些。比如由圖可見第5列的第6個矩陣元的數值大于第4個矩陣元，因此C5→C6的電子轉移量大于C5→C4的。
我們再看看其它激發態的熱圖。S0→S2的第1種模式構建的原子ATM圖，以及S0→S9的原子-原子電荷轉移矩陣圖如下所示，對應的空穴&電子等值面圖也都附上了。由于S0→S9明顯牽扯到了氫原子，因此繪圖的時候也把氫考慮了進去（之所以S0→S9激發沒有用原子TDM圖展示，是因為此圖與空穴&電子等值面圖對應并不很理想）
S0→S2的熱圖右上角有一塊數值較大區域，體現出在體系末端的硝基部分，空穴和電子同時有較大分布。而且圖像最右側一列數值都不是很小，因此可以認為硝基向體系的中央區域轉移了一定量的電子，這與空穴&電子等值面圖上看到的情況相一致，也可以描述為硝基部分與體系中間區域在S0→S2激發時存在所謂的“相干”。
從S0→S9的原子-原子電荷躍遷矩陣熱圖中可以看出，電子從1~5、7~9號原子的區域往13號氫原子上轉移電子非常明顯，這和空穴&電子等值面圖傳遞的信息完全對應。6號原子上幾乎只有明顯的綠色等值面，說明6號原子的電子沒有給其它原子，而是從其它原子上獲得了不少，相應地，熱圖上Y=6那一行顏色也比較鮮明，而X=6那一列則顏色很暗。可見，空穴&電子等值面圖視覺效果最為直觀，但是如果結合原子-原子電荷躍遷矩陣熱圖一起討論，可以從定量角度理解得更透徹，而且也避免了等值面數值設定的任意性導致做出不合理判斷的可能。
3.3 繪制基態到激發態的片段TDM的熱圖
這一節我們用下面這個體系作為例子，主要演示一下繪制基于片段的TDM的熱圖。用到的文件可以在這里下載：http://www.shanxitv.org/attach/436/file.rar。
繪制片段TDM的熱圖的時候片段可以隨意定義，矩陣元的序號就是片段的序號，所以也就不用要求體系中的原子序號必須和體系中原子坐標有對應關系，因此對于任意形狀的體系，哪怕是星形、環形等都可以使用，非常靈活，而且圖像展現的信息比原子TDM熱圖更為緊湊。當前我們按照上圖標注的顏色對片段進行劃分。
此體系的Gaussian的TDDFT輸入文件是文件包里的tdmat.gjf，關鍵詞為# CAM-B3LYP/6-31G* TD(nstates=10) IOp(9/40=4)。運行后得到文件包里的tdmat.out，同時產生的chk文件轉換后是文件包里的tdmat.fchk。
繪制片段TDM熱圖前我們先看看S0→S1的原子TDM熱圖。啟動Multiwfn，依次輸入
tdmat.fchk
18 //電子激發分析
2 //繪制原子/片段躍遷矩陣的熱圖
tdmat.out
1 //考察S0→S1激發。然后程序會產生這個躍遷對應的TDM
n //不對剛產生的TDM做對稱化
1 //用第1種模式將TDM收縮為原子TDM
1 //繪制原子TDM的熱圖
看到的圖像如下
像當前體系這樣原子數比較多，又含有環的情況，將體系結構和熱圖上的坐標進行對應比較麻煩，雖然也可以自行編輯圖片來標記坐標軸什么范圍對應體系的什么區域，比如像我很早前寫的http://www.shanxitv.org/136文中的圖那樣，但比較費事，而且如果體系特征更復雜的話，很難去標記。若是原子序號都不是連貫的，甚至根本沒法標記。而自定義片段來繪制片段TDM就沒這個問題了。
片段可以在Multiwfn界面里直接輸入，也可以從文本文件中讀取。如果片段多的話，建議還是寫文本文件來定義片段，這樣以后重新作圖時就不用再次重新輸入一遍了，也免得輸入的時候萬一輸錯而前功盡棄。記錄片段設定的文件名隨意，我們新建一個叫tdmfrag.txt的文件，寫上以下內容，每一行定義一個片段：
1-23
24-33
34-43
44-55
56-63
橫杠代表范圍，格式很靈活，也可以寫成例如2,4,9-14,88這種更復雜的形式。體系比較大的時候，如果不好確定片段里的原子序號，可以用gview來輔助給出，做法參看本文文首的視頻里的演示。
然后在Multiwfn窗口里繼續輸入
-1 //讀入片段定義
0 //從外部文件中讀入
tdmfrag.txt //定義了片段的文件路徑
5 //修改色彩刻度范圍
0,0.4 //把色彩刻度下限和上限設為0和0.4
1 //繪制熱圖
此時看到下圖，橫/縱坐標序號現在對應的是片段序號了，空穴&電子等值面圖也一并附上（isovalue用的是0.001），藍框標注的部分是4號片段（己三烯）
根據熱圖中的顏色，我們得知電子和空穴大部分都在片段4上，也一定程度同時出現在片段1和5上，這和等值面圖看到的情況大體一致。由于圖上非對角元不大，因此這種電子激發沒有造成很顯著的片段間電子轉移。粗略地說，這個激發的主體特征是片段4上的局域激發。
如果進入熱圖繪制功能的時候載入的是IFCT模塊產生的atmCTmat.txt文件，然后再同上載入片段設定，那么繪制的就是片段-片段電荷轉移矩陣的熱圖了。這里就不再演示了，請大家自行繪制。
下面把S2~7的TDM熱圖繪制出來并手工合并在一起，得到下圖：
S2~S5的圖的非對角元的數值相對于對角元來說都不顯著，這種情況下只要看哪些對角元較大，就知道哪些片段在電子激發過程中被明顯牽扯到。可見S0到S2、S3的躍遷主要都是在片段1，即分子中蒽的那部分，但S2也稍微牽扯到片段4，總的來說這兩種躍遷都可以算是局部激發。而S0→S4一看就是整體激發。S5可以判斷為主要是片段3那個苯環的局域激發，但也多多少少牽扯到其近鄰區域。S6和S7彼此間有點呈現鏡像關系，由圖可見體系各個片段在電子激發時都牽扯到了，要么有空穴分布，要么有電子分布，要么二者都有。對于S0→S6來說，根據熱圖的顏色我們可以推測片段2、3、4向片段1轉移了電子，片段3、5也往片段4有電子轉移。如果你寫文章的時候要同時討論基態到一大批激發態的特征，顯然像這樣把相應的熱圖都合并到一起放到文章里是非常直觀的。
3.4 繪制躍遷偶極矩矩陣的熱圖
這一節演示一下怎么繪制躍遷偶極矩矩陣的熱圖，考察這樣的圖有助于我們搞清楚體系不同區域對躍遷偶極矩的貢獻，它和躍遷偶極矩密度等值面圖有類似的價值，但表現形式不同。這里還用NH2-C8-NO2那個體系的S0→S1躍遷作為例子。
首先生成原子的躍遷偶極矩矩陣。啟動Multiwfn，依次輸入
examples\excit\NH2_C8_NO2\NH2_C8_NO2.fchk
18 //電子激發分析
11 //將躍遷電/磁偶極矩分解為基函數和原子的貢獻
examples\excit\NH2_C8_NO2\NH2_C8_NO2.out
1 //考察S0→S1激發
1 //電偶極矩
y //將原子的躍遷偶極矩矩陣導出
如屏幕的提示所示，矩陣數據已經導出到當前目錄下AAtrdip開頭的文件中了。其中AAtrdipX.txt是原子躍遷偶極矩矩陣X分量的數據文件，本例我們就是要考察躍遷偶極矩的X分量的情況，因此下面要用這個文件。
重啟Multiwfn，輸入
o //載入上次的輸入文件
18 //電子激發分析
2 //繪制原子或片段的躍遷矩陣的熱圖
AAtrdipX.txt
此時屏幕上看到
Sum of all elements (including hydrogens): -4.38875223
Maximum and minimum (including hydrogens): 0.65818572 -0.82543408
Sum of all elements (without hydrogens): -2.86353889
Maximum and minimum (without hydrogens): 0.65818572 -0.82543408
此處的-4.38875223是包括氫在內的所有矩陣元的加和。由于躍遷偶極矩矩陣某個分量的所有矩陣元加和正對應于躍遷偶極矩的這個分量，因此-4.38875223就是S0→S1激發的躍遷偶極的X分量，這和Gaussian輸出文件里下面這部分看到的-4.3881是直接對應的。
Ground to excited state transition electric dipole moments (Au):
state X Y Z Dip. S. Osc.
1 -4.3881 -0.2570 0.0075 19.3211 1.6097
如上面提示所示，當前的矩陣元最小值是個負值-0.825，但是繪制熱圖的功能默認把色彩刻度下限設為0，因此我們得改一下色彩刻度，而且最好要讓色彩刻度上、下限的絕對值相同，這樣正負部分刻度是對稱的。可以反復嘗試找到令圖像最能充分反映矩陣特征的數值。范圍設太小的話，超過色彩刻度上限和下限的部分分別會顯示為白色和黑色，就不美觀了；而范圍太大的話，矩陣元的差異就不容易通過顏色分辨開。
我們接著在Multiwfn里輸入
5 //修改色彩刻度范圍
-0.7,0.7 //色彩刻度下限和上限
1 //繪制熱圖
所得圖像如下。這里也把這個激發的躍遷偶極矩密度X分量的等值面圖一起附上便于對比，等值面數值設為了0.01
這個熱圖越藍(越紅)的矩陣元對躍遷偶極矩X分量起到越負(越正)的貢獻。由于圖上大部分都是藍色的，因此所有矩陣元加和必然為負，解釋了躍遷偶極矩X分量為什么為顯著的負值-4.388。圖上明顯偏離對角線的矩陣元數值都很接近0，即都是綠色，因此原子間的長程耦合對于躍遷偶極矩X分量基本沒有貢獻。圖上有幾個部分藍色很顯著，比如(2,2)附近、(9,9)附近等，體現出這些位點及近鄰的區域產生了明顯的負貢獻，這也正好對應等值面圖上這些區域以藍色等值面為主的現象。也有的位點比如(1,2)是明顯的正值，這說明這倆原子間的耦合對躍遷偶極矩X分量有明顯正貢獻，這也正對應于等值面圖上所看到的1-2原子間都是綠色等值面的現象。熱圖的中部基本都是一片綠，數值很小，相應地等值面圖上分子中部沒有什么等值面出現。
可見，將躍遷偶極矩密度和躍遷偶極矩矩陣結合到一起分析，對于弄清楚躍遷偶極矩的內在構成很有幫助。
3.5 繪制激發態間的躍遷密度矩陣的熱圖
這一節演示用Multiwfn繪制激發態到激發態的躍遷密度矩陣的熱圖，用NH2-C8-NO2的S1→S2躍遷作為例子。
首先產生S1→S2對應的躍遷密度矩陣數據文件。啟動Multiwfn，依次輸入
examples\excit\NH2_C8_NO2\NH2_C8_NO2.fchk
18 //電子激發分析
9 //產生躍遷密度矩陣
2 //產生兩個激發態間的躍遷密度矩陣
examples\excit\NH2_C8_NO2\NH2_C8_NO2.out
1,2 //考察S1→S2激發
直接按回車用默認的閾值
n //不對產生的躍遷密度矩陣做對稱化
n //不產生TDM.fch文件
現在當前目錄下已經有了記錄S1→S2躍遷密度矩陣的tdmat.txt。
重新啟動Multiwfn，輸入
o //載入上次的輸入文件
18 //電子激發分析
2 //繪制原子或片段的躍遷矩陣的熱圖
tdmat.txt
1 //以方式1構建原子躍遷密度矩陣
1 //繪制熱圖
得到的圖像如下。S2的密度減去S1的密度對應的密度差的等值面圖也一并附上了，怎么繪制看《使用Multiwfn計算激發態之間的密度差》（http://www.shanxitv.org/429）。
由于熱圖右側X=11,12那一條中1~9號原子對應的區域數值很大，因此我們可推測在S1→S2激發過程中11、12號原子向1~9號原子轉移了大量電子，這和等值面圖上看到的結論完全一致。此例體現出密度矩陣熱圖不僅對于分析基態到激發態的躍遷很有益，對于討論假想的激發態之間的躍遷也同樣可以用。
4 技巧：對一批激發態批量繪制熱圖
如果你一次要考察一大批激發態，一個一個生成圖像文件嫌麻煩，那么可以用Linux的shell腳本實現。
一次性產生指定范圍激發態的片段TDM熱圖的腳本是Multiwfn程序包里examples\scripts\allTDM.sh文件。把3.3節用到的tdmat.fchk、tdmat.out、tdmfrag.txt和allTDM.sh都放到Multiwfn目錄中，然后進入Multiwfn目錄，用chmod +x allTDM.sh給其加上可執行權限，然后運行./allTDM.sh，這個腳本就會自動調用Multiwfn，在當前目錄產生1.png、2.png直到7.png，對應S0-S1、S0-S2直到S0-S7的以模式1方式構建的片段TDM的熱圖。整個運行過程轉眼就能完成，極其方便。
一次性產生指定范圍激發態的片段-片段電荷轉移矩陣的熱圖的腳本是Multiwfn程序包里examples\scripts\allCTmat.sh，執行方式同上。
默認的輸入文件名、要算的激發態范圍都可以自己編輯腳本來修改，腳本很簡單，一看就秒懂。如果對如何通過腳本批量執行Multiwfn一無所知的話，建議看手冊5.3節的相關說明。
5 總結&其它
本文全面地介紹了躍遷密度、躍遷密度矩陣、躍遷偶極矩密度、躍遷偶極矩矩陣、電荷轉移矩陣等概念，并結合實例介紹了如何通過Multiwfn進行繪制和分析，同時還與Multiwfn的空穴-電子分析、密度差圖等進行了對照。例子充分展現出分析上述函數和矩陣對于表征電子激發、弄清楚電子激發以及躍遷偶極矩的內在特征非常有幫助，十分值得在涉及電子激發問題的研究中充分靈活運用。
最后再強調一下，TDM和電荷轉移矩陣雖然大多數情況下定性一致，但是有時也會出現定性不同的情況，尤其是非對角元上的差異可能較大，比如3.3節例子里的S0→S4激發。存在明顯差異的時候，一般都是電荷轉移矩陣元數值較大而TDM矩陣元較小（更深層的原因一定程度上是因為計算TDM非對角元的時候不同組態對應的項出現一定程度的正負抵消）。當電子激發沒有一個組態函數占明顯主導的情況，這兩種矩陣出現定性不符的幾率相對較大。凡是遇到把空穴&電子等值面圖和TDM放在一起時發現特征不怎么對應的時候，建議使用電荷轉移矩陣代替TDM。鑒于電荷轉移矩陣物理意義更清楚，只用它來討論而不用TDM是完全可以的。
有人問對比一系列體系或者一個體系的不同激發態的時候是否需要把TDM、電荷轉移矩陣的色彩刻度統一。在我來看，如果你的目的是判斷各個激發態的特征，比如是LE還是CT、涉及哪些片段間電荷轉移，那么色彩刻度不是必須統一，讓各個態的特征都能展現得盡可能鮮明才好。而如果你要定量考察不同躍遷中片段間耦合、電子轉移的程度差異，那么色彩刻度應當統一。
在《一篇最全面、系統的研究新穎獨特的18碳環的理論文章》（http://www.shanxitv.org/524）中介紹的筆者的論文中（后來其中電子激發和非線性光學部分專門正式發表于Carbon, 165, 461-467 (2020)），筆者使用躍遷偶極矩密度方法探究了18碳環這種電子結構十分特殊的體系的S0->S21激發具有非常強吸收的原因，是躍遷偶極矩密度分析的很好的范例，非常推薦讀者一讀。
本文中各種躍遷矩陣圖的顏色都是以默認的彩虹方式變化的，色彩變化實際上是可以在Multiwfn的繪制躍遷矩陣的功能里通過選項“9 Set color transition”隨意調節的，支持一大堆色彩變化方式，比如藍-白-紅、灰度等等。Multiwfn產生的矩陣數據也可以導出，然后放到第三方作圖程序如Origin里繪制成熱圖，屆時有更多選項可以調節作圖效果，更為靈活。在Origin里的繪制方法見本文頂端的視頻里的演示。