頻域空間密度矩陣重正化群的研究進展

姜 童，任佳駿，帥志剛

（清華大學(xué)化學(xué)系,教育部有機光電子與分子工程重點實驗室,北京100084）

在量子多體問題中，波函數(shù)的精確求解受到“指數(shù)墻”的制約，因此在有限的計算資源下對波函數(shù)進行合理的近似是一種常用手段.1992年，White［1，2］提出了密度矩陣重正化群（Density matrix renormalization group，DMRG）理論.DMRG利用約化密度矩陣的本征矢量在實空間的局域性來降低希爾伯特空間的維數(shù)，通過反復(fù)迭代變分優(yōu)化系數(shù)，成為了一種求解多體波函數(shù)的有效方法.該方法成功地應(yīng)用于求解一維強關(guān)聯(lián)格點模型中的低能電子態(tài)的電子結(jié)構(gòu)，如對于短程作用的Hubbard模型、Heisenberg模型等給出了接近精確解的結(jié)果.Shuai等［3～8］很快地將DMRG推廣到具有長程相互作用的半經(jīng)驗量子化學(xué)中，并用以研究共軛高分子的低能級激發(fā)態(tài)次序以及線性、非線性光學(xué)響應(yīng)等問題，這被國際學(xué)術(shù)界視為DMRG在量子化學(xué)領(lǐng)域應(yīng)用的開端.White等［9］和Chan等［10］發(fā)展了從頭算的DMRG算法，并成功應(yīng)用于過渡金屬中心配合物以及大有機共軛體系的電子結(jié)構(gòu)計算中.DMRG的含時理論和響應(yīng)理論（Time dependent DMRG，TD-DMRG）［11～16］以及頻率空間算法［5，8，17～27］的快速發(fā)展使得DMRG能夠用于動態(tài)響應(yīng)性質(zhì)的計算.TD-DMRG通過含時演化得到時間關(guān)聯(lián)函數(shù)，而后通過傅里葉變換得到頻域信息，該方法在包括線性光譜計算［15］、激子分離［28］、分子的超快內(nèi)轉(zhuǎn)換［29，30］以及單態(tài)裂分［31，32］等有機半導(dǎo)體體系的光物理過程得到了廣泛應(yīng)用.然而在含時演化過程中體系的糾纏會隨著時間快速增長［33～35］，在有限的變分空間中演化的誤差增大，從而對長時間演化帶來了挑戰(zhàn)，另一方面，傅里葉變換得到的譜圖的精度取決于演化時間的長度.

頻率空間的算法為動態(tài)響應(yīng)函數(shù)的計算提供了不同角度.1995年，Hallberg［17］首次將基于Lanczos遞歸的連分?jǐn)?shù)展開（Continue fraction expansion，CFE）方法與DMRG結(jié)合，并應(yīng)用于頻域自旋關(guān)聯(lián)函數(shù)的求解.這種方法簡單易行，但精度較低，只適用于計算低能態(tài)的較為簡單的分立譜［20］.隨后，Ramasesha和Shuai等［5，8］提出了更精確的修正矢量密度矩陣重正化群CV-DMRG（Correction vector DMRG）方法，可以精確地計算復(fù)雜的連續(xù)譜函數(shù)，并適用于非線性光學(xué)響應(yīng)函數(shù)，但需要的計算量更大.基于CV-DMRG的工作方程，Jeckelmann［21］提出一種變分優(yōu)化的動態(tài)DMRG算法（Dynamical DMRG），可以將計算精度提高.該方法更為簡潔有效，被視為最精確的計算譜函數(shù)的DMRG算法［14，25］.最近，Schollw?ck和Delft等［23］將矩陣乘積態(tài)（MPS）與Chebyshev多項式展開相結(jié)合提出了CheMPS算法，該方法在自旋模型［23］及雜質(zhì)模型［36］的譜函數(shù)計算中實現(xiàn)了計算量與精度之間較好的平衡，但其潛力仍有待進一步探索.Chan等［25～27］和Haegeman等［37］基于含時微擾理論和DMRG波函數(shù)的切空間，提出了解析線性響應(yīng)DMRG方法，也被用來求解線性響應(yīng)性質(zhì).本文將介紹這些頻率域發(fā)展的動態(tài)DMRG響應(yīng)理論，并重點介紹我們課題組在有限溫度發(fā)展的算法.

1 密度矩陣重正化群與矩陣乘積態(tài)

對于含有N個軌道的系統(tǒng)（i=1，…N），|σi〉代表第i個軌道的d維的正交基組（如對于電子系統(tǒng)，每個空間軌道存在|↑〉，|↓〉，|↑↓〉，|vacuum〉4種電子占據(jù)狀態(tài)，d=4；而對于玻色子系統(tǒng)，d原則上是無限大，因此需要做截斷），整個系統(tǒng)的希爾伯特空間|σ1…σN〉=|σ1〉?…?|σN〉為dN維，精確解可以用完全組態(tài)相互作用（Full configuration interaction，F(xiàn)CI）的波函數(shù)表示：

式中，系數(shù)cσ1…σN共有dN種取值，因此對于稍大的體系便難以實現(xiàn)精確對角化.DMRG波函數(shù)采用可變分優(yōu)化的MPS作為波函數(shù)的擬設(shè)，即將高維向量cσ1…σN分解為N個低秩矩陣乘積的形式［38，39］：

式中，Lσ與Rσ分別是左規(guī)整與右規(guī)整的矩陣，即［見圖1（D）和（E）］為重正化基|lai-1〉|σi〉|rai〉下的系數(shù)矩陣：

Fig.1 Graphical representation of matrix product state/operator(MPS/MPO)of wavefunction and operator

傳統(tǒng)的DMRG并不顯式地構(gòu)造全空間下的算符，而是處理算符在重正化基下的投影，如第i個格點的哈密頓量的投影其中在MPS出現(xiàn)后，對于任意的算符也可以被自然地表示為矩陣乘積算符（Matrix product operator，MPO）的形式（如下式），見圖1（B），MPO的局域矩陣是一個四維張量，wi-1與wi被稱為“虛擬鍵”被稱為“物理鍵”.這使得DMRG在形式上變得更為簡潔，同時也提高了某些算法的精度.

Fig.2 Adding two different MPOs

當(dāng)一個局域矩陣維數(shù)為D的MPO作用到虛擬鍵維數(shù)為M的MPS上后，會得到一個維數(shù)為MD的新MPS，如圖3所示，在很多場景中，新得到的MPS的維數(shù)較高，因此需要通過奇異值分解（SVD）分解或者變分的方式做壓縮.

Fig.3 Applying an MPO to an MPS

Fig.4 Overlaps and expectation values

DMRG通過掃描的方式對波函數(shù)進行迭代優(yōu)化，以單格點的DMRG算法為例，當(dāng)優(yōu)化到第i個格點，對式（3）中的Cσi變分來最小化拉格朗日函數(shù)為變分系數(shù)，其收斂的結(jié)果對應(yīng)于的本征能量，而后得到式（7），利用的左右規(guī)整性，可將該方程進一步簡化為一個特征值問題，見圖5.

Fig.5 Eigenvalue problem when minimizing the Lagrangian with respect to local matrix Cσi

2 線性響應(yīng)理論

下面將簡要介紹線性響應(yīng)理論，并推導(dǎo)頻率空間的響應(yīng)函數(shù)在有限溫度及零溫度下的函數(shù)形式.未加說明的，物理量均用原子單位.

對于處于熱平衡（溫度為T）的正則系綜，其密度矩陣表示為其中配分函數(shù)Z(β)=為玻爾茲曼常數(shù)當(dāng)系統(tǒng)受到來自外界勢場的f(t)V的擾動時，其哈密頓量表示為

式中，f(t)為外勢場的強度；V^為系統(tǒng)的可觀測物理量，滿足V^=V^?；當(dāng)系統(tǒng)受到的擾動比較小時，利用線性響應(yīng)理論，O^(t)的熱力學(xué)期望值用Kubo方程［40］表示為

式中，χ(t-t′)被稱為推遲格林函數(shù)（僅當(dāng)t＞t′時，不為0），

式中，θ(t-t′)為單位階躍函數(shù)；O^(τ)=eiH^0τO^e-iH^0τ；[·，·]為對易關(guān)系（即是指關(guān)于的密度矩陣求期望.將χ(t-t′)做傅里葉變換到頻率空間可得：

式中，Em，En分別對應(yīng)于的第m和第n個本征態(tài)的本征能量.

式中，

0 為基態(tài)波函數(shù).

其正比于體系吸收外界勢場能量的速率.

3 零溫度下的算法

3.1 Lanczos DMRG

Lanczos迭代算法在數(shù)學(xué)上被用來求解大型稀疏厄米矩陣的線性方程與特征值問題，其思想是將原矩陣投影到Lanczos遞歸產(chǎn)生的向量所張成的Krylov子空間上，從而大大降低對角化的難度.該方法最早被Gagliano與Balseiro［42］用來求解零溫度下的基態(tài)與線性響應(yīng)，Hallberg［17］將其與DMRG相結(jié)合來計算T=0的單粒子響應(yīng)函數(shù)（下式），使其能夠處理更大的系統(tǒng)［18，19，43～46］.

式中，第二個等號后為萊曼表示的譜函數(shù).Lanczos DMRG將難以直接處理的哈密頓量投影到Krylov子空間的三對角矩陣Heff：

直接對角化得到子空間中的特征向量與特征值，代入式（15）即可.Heff由Lanczos迭代算法得到，

除了直接對角化Heff，通過連分?jǐn)?shù)展開也可直接得到Lorentz展寬為η的光滑譜圖：

式中，z=E0+ω+iη.

Lanczos算法可以在DMRG基態(tài)程序基礎(chǔ)上較為容易地實現(xiàn)，同時計算量小，在MPS/MPO的框架下，Lanczos算法只需進行將MPO作用到以及MPS之間的求和操作即可快速實現(xiàn)［19］，其實現(xiàn)的細(xì)節(jié)請參見第1節(jié).但Lanczos只能應(yīng)用于計算低能態(tài)的較為簡單的分立譜［20］，這是由于在迭代過程中存在誤差累積，同時由于第n+1個Lanczos向量僅依賴于第n個和第n-1個向量，因此與之前n-2個向量之間的正交性逐漸被破壞，Kühner和White［20］提出作為判據(jù)及時終止迭代.Dargel等［19］提出了將Lanczos向量重新做正交化的方法，以改善該問題.

3.2 Correction Vector DMRG

早在DMRG被提出之前，Soos與Ramasesha［47］就提出了在精確對角化框架下的修正矢量（Correction vector，CV）方法，求解π電子共軛體系的Pariser-Parr-Pople（PPP）模型的線性與非線性響應(yīng)光學(xué)性質(zhì).該方法盡管精確，但限于較小的體系（N＜12）.之后，Ramasesha等［5，8］將DMRG與CV相結(jié)合，能夠處理更大的體系.為求解式（15），將頻率為ω時的CV定義如下：

式（22）可通過掃描的方式在局部的重正化基上利用共軛梯度法進行求解，當(dāng)采用的展寬η較小時，線性方程在達(dá)到共振的位置是接近奇異的，為了提高收斂速度，可以采用預(yù)處理條件來減小矩陣的病態(tài)程度.在重正化基表示中，第i個格點的投影哈密頓量為需要說明的是，式（22）中的引入增大了條件數(shù)（約為原條件數(shù)的平方），增加了系統(tǒng)的病態(tài)程度，給求解帶來困難.此外，傳統(tǒng)的DMRG算法不能直接給出需要通過對其進行近似，二者在M→+∞時嚴(yán)格相等，但是在有限的M時會引入額外的誤差.

這種方法避免了直接求解式（22），與原始Lanczos DMRG不同的是，這里的密度矩陣考慮了的貢獻，即對于不同的ω采用了不同的重正化基去計算，因此更為準(zhǔn)確.

3.3 Dynamical DMRG

Jeckelmann［21］基于CV-DMRG進一步提出了變分優(yōu)化的方法，稱為動態(tài)DMRG（Dynamical DMRG，DDMRG），將求解CV-DMRG的線性方程［式（22）］轉(zhuǎn)化為最小化的問題：

當(dāng)F處于極小點時，對于譜函數(shù)S(ω)［見式（21）］，有因此不必顯式地利用求解S(ω)，只要求得F的最小值即可.在此框架下，對于給定的若其誤差為O(ε)，則S(ω)的計算誤差變?yōu)镺(ε2)，這便是最小化泛函方法的優(yōu)勢所在.在掃描優(yōu)化F的過程中，對第i個格點矩陣Aσi進行變分，其余格點保持固定，此時的工作方程為

DDMRG的另一項優(yōu)勢在于其可以自然地推廣到MPS/MPO的框架下，式（26）在MPS/MPO框架下的表示如圖6（A）所示.這一線性方程可以通過標(biāo)準(zhǔn)的迭代算法，如共軛梯度法進行求解.圖6（B）表示將求解后的解進行奇異值分解后，更新局域矩陣，并為下一個位點提供初猜.在傳統(tǒng)的DDMRG算法中與基態(tài)波函數(shù)通過態(tài)平均共用一套重正化基［21］，通過將用2個不同的MPS分別表示，將每個態(tài)都表示得更加準(zhǔn)確，從而提高了計算精度.另一方面，在以往的算法中，局部的表示在重正化基上（見3.2節(jié)），因此需要對于包括在內(nèi)的高階矩進行近似，而利用MPO表示則可以對此類高階矩進行精確表示.

Fig.6 Graphical representation of minimizing the Lagrangian with respect to local matrix Aσi[22]

3.4 Chebyshev多項式展開矩陣乘積態(tài)

類似于Lanczos方法，Chebyshev多項式展開的方法利用遞歸產(chǎn)生一系列的Chebyshev向量來展開響應(yīng)函數(shù)，并被成功應(yīng)用到Heisenberg模型、Holstein模型、Anderson模型等在零溫度或有限溫度下的光吸收譜、光電導(dǎo)率等響應(yīng)性質(zhì)的計算［49～52］.Schollw?ck和Delft等［23］提出了將Chebyshev多項式展開與MPS結(jié)合的CheMPS方法.Chebyshev展開的核心思想是采用核多項式展開去近似線性響應(yīng)函數(shù)中的δ函數(shù)：

首先，介紹一下其數(shù)學(xué)背景［52］，對于一般的函數(shù)f(x)，可以通過Chebyshev多項式進行展開：

其中，Tn(x)服從三項遞推關(guān)系：

上述展開僅在x∈[-1，1]成立.因此，為了處理首先做投影與在實際操作中，為了避免數(shù)值不穩(wěn)定帶來的問題，要求ω′∈[-W′，W′]，其中W′略小于1：

式中，gn是阻尼系數(shù)，這是由于在實際操作中多項式不能展開到無窮階，因此會產(chǎn)生所謂的吉布斯振蕩［23］，引入gn可以達(dá)到平滑光譜的效果，如最常用的Jackson系數(shù)與Lorentz系數(shù)可分別在頻率空間達(dá)

將式（31）代入式（27）可得：

由于系統(tǒng)的整個能量區(qū)間W可能很大，此時為了得到較高精細(xì)程度的光譜需要演化很多的步數(shù)，而實際計算上，存在響應(yīng)的頻率區(qū)間相對于系統(tǒng)的整個能量空間會較小，因此可選擇將一個有效的能量區(qū)間W*投影到[-1，1]之間.當(dāng)W*＜W時，在每次演化得到對其進行壓縮的過程中會引入“高能態(tài)”（對應(yīng)于的特征值大于1的分量），這將導(dǎo)致接下來的演化發(fā)散，在此情況下，就很有必要進行額外的能量截斷過程［23］，需要通過能量截斷將每次引入的“高能態(tài)”扔掉.由于的希爾伯特空間難以直接處理，因此需要在局域的格點進行操作，在第i個格點對應(yīng)的重正化基下的下構(gòu)造維度為dk的Krylov子空間，將在子空間下的矩陣做對角化可以選取出能量的絕對值大于1的分量構(gòu)造投影算符而后作用到的第i個格點，移動到下一個點，掃描ns圈后結(jié)束能量截斷.需要注意的是，并沒有嚴(yán)格的方法可以將高能部分全部扔掉，ns與dk的最優(yōu)組合也并不直接可以得到，ns與dk需要經(jīng)驗地進行選擇.

3.5 解析的線性響應(yīng)DMRG

上述算法可歸為一類，首先，這幾種方法都直接面向求解線性響應(yīng)函數(shù)［見式（15）］，其次，在與DMRG結(jié)合之前，Lanczos方法、CV方法以及Chebyshev多項式展開方法便被獨立地用來直接求解線性響應(yīng)函數(shù).與上述方法不同，解析的線性響應(yīng)DMRG方法基于DMRG的矩陣乘積態(tài)結(jié)構(gòu)，通過微擾理論直接計算波函數(shù)的局域格點矩陣對于外界勢場微擾的響應(yīng)［25～27，37］.

對于第1節(jié)中關(guān)于系統(tǒng)波函數(shù)的混合規(guī)整表示，將基態(tài)波函數(shù)表示為

式中，上角標(biāo)(0)代表系統(tǒng)未經(jīng)外界勢場擾動.當(dāng)系統(tǒng)受到外界光場V(t)=Veiωt+V*e-iωt的微擾，使得各格點矩陣產(chǎn)生響應(yīng)：

接下來，移動到第i+1個格點求解當(dāng)?shù)玫绞諗康囊浑A線性響應(yīng)波函數(shù)后，可以求解線性極化率（ω＞0）：

G(ω)對應(yīng)于式（13）中的第一項，其中μ(0)i是偶極算符在第i個格點對應(yīng)的重正化基下的投影，C(1)(ω)是關(guān)于頻率為ω的光場擾動的一階響應(yīng)波函數(shù).

4 有限溫度下的算法

4.1 有限溫度密度算符的求解

DMRG本身是針對波函數(shù)的理論，對于有限溫度或混合態(tài)的情況下則需要用密度矩陣（算符）描述，在熱場動力學(xué)算法（Thermal field dynamics，TFD）中［53］，處于溫度T(T＞0)的任意量子態(tài)可以表示在由物理空間P與輔助空間Q構(gòu)成的擴大的希爾伯特空間P?Q中（此輔助空間Q是物理空間P的一個拷貝）.

Fig.7 Graphical representation of obtaining the ensemble expectation value of operator(marked in blue)by thermal field dynamics

因此，ρβ可以表示為在DMRG的框架下，熱場動力學(xué)也被稱為純化的方法［38，54～56］，借助于MPS/MPO的語言，可以對這一過程進行更為直觀的理解，由式（44），可以在虛時間軸從τ=0演化到τ=β/2得到：

由于式（46）不是幺正的，在每演化一步后，ρ(τ)通過進行歸一化.

4.2 Lanczos DMRG

Prelov?ek［57］等從另一個角度通過對初態(tài)進行采樣的方式來近似得到熱平衡密度矩陣，從而將Lanczos DMRG擴展到了有限溫度.T＞0時的密度矩陣表示為

式中，R為取樣的數(shù)目.對做Lanczos迭代得到三對角的Heff，對角化得到一組特征向量以及特征值從而將作用于上近似為

對于有限溫度下的自相關(guān)函數(shù)χ(1)(ω)［見式（12）］：

4.3 Dynamical DMRG

最近，Shuai等［22］在MPS的框架下，通過純化密度矩陣的方法建立了有限溫度下的DDMRG.對于有限溫度下的響應(yīng)函數(shù)S(ω)（參照第2節(jié)）：

基于此構(gòu)造泛函F：

當(dāng)X(ω)滿足式（55）時，F(xiàn)取得全局最小值，同時滿足S(ω)=-Fmin/πη.X(ω)可以通過迭代求解方程（下式）進行變分優(yōu)化［38］（如圖8所示）：

Fig.8 Linear equation?F/?Aσi=0（Eq.55）to optimize the local site Aσ（iin purple）at finite temperature［22］is marked in red.The operators are shown in blue.Copyright 2020，American Chemical Society.

4.4 Chebyshev展開矩陣乘積態(tài)

在3.4節(jié)介紹了零溫度下采用Chebyshev展開的DMRG方法.Tiegel等［24］將其推廣到求解有限溫度的譜函數(shù)：

如同零溫度的做法，需要先將ω與L投影，對于ω∈[ωmin，ωmax]，作如下操作：

其中，W=ωmax-ωmin，

將式（60）代入式（58），可得到與在零溫度下相同的表達(dá)式［參見3.4節(jié)中式（34）］，其中Chebyshev矩滿足三項遞推關(guān)系：

5 應(yīng) 用

自1995年Hallberg［17］提出Lanczos DMRG求解一維Heisenberg模型的自旋關(guān)聯(lián)函數(shù)以來，頻域空間的DMRG算法逐漸發(fā)展，并被廣泛應(yīng)用于電子體系（Hubbard模型及其擴展、Anderson雜質(zhì)模型、從頭算量子化學(xué)哈密頓量）、電-聲子體系（Holstein模型）、自旋體系（Heisenberg模型）等不同體系響應(yīng)性質(zhì)的計算.下面將從電子問題以及電聲子耦合問題兩個方面介紹其代表性的應(yīng)用場景.

5.1 電子問題

無機半導(dǎo)體一般具有較大的介電常數(shù)ε＞10，因此對于固相下的分子間庫侖相互作用可以實現(xiàn)很好的屏蔽，故常用單電子的能帶模型.然而，對于有機半導(dǎo)體而言，分子間的相互作用力以相對小的范德華力為主，一般具有較小的介電常數(shù)（ε約為2～4），對于分子間的庫侖相互作用不能較好地屏蔽，體系經(jīng)光激發(fā)產(chǎn)生的電子-空穴對（激子）的結(jié)合能較大，激子相對局域，半徑較小.擴展的Hubbard-Peierls模型是考慮電子-空穴之間的吸引勢的一種簡單處理方式：

式中，t為最近鄰跳躍積分；δ為調(diào)節(jié)鍵長的參數(shù)，描述了Peierls不穩(wěn)定性.如對于聚乙炔鏈來說，雙鍵對應(yīng)的跳躍積分為t（1+δ），單鍵對應(yīng)的跳躍積分為t（1-δ），U為占位庫侖排斥，V為最近鄰的庫侖勢，如圖9所示.

基于此哈密頓量，Pati等［8］研究了不同長度的鏈在不同參數(shù)下的三階非線性極化率，見圖10.3.2節(jié)已經(jīng)介紹了基于CV-DMRG如何求解線性響應(yīng)性質(zhì)，盡管非線性響應(yīng)的表示更為復(fù)雜，但仍可采用CV-DMRG在同一框架下進行計算［5，8］.

Fig.9 Representation of extended Hubbard-Peierls model

Fig.10 Plot of the log of average third-order polarizability log of the chain length L[8]

Jeckelmann［21］基于CV-DMRG方法提出了DDMRG方法，并通過計算流-流關(guān)聯(lián)函數(shù)計算上述模型在不考慮近鄰電子關(guān)聯(lián)（V=0）時的光電導(dǎo)譜（見圖11），當(dāng)再忽略式（62）中的Hubbard項（U=0）時，此哈密頓量描述了自由電子體系，因此其光電導(dǎo)率可以精確求解，由圖11（A）可見，基于DDMRG的結(jié)果與精確解完全重合，印證了該方法的精確性和有效性.圖11（B）描述了對于Mott-Hubbard模型、Peierls模型以及Hubbard-Peierls模型3種不同參數(shù)條件下的光電導(dǎo)率.

Fig.11 Optical conductivity calculated by DDMRG［21］

Chan等［25］發(fā)展了解析的線性響應(yīng)DMRG方法，通過計算聚乙炔的極化率，將其與DDMRG方法進行了比較.對于較小的M（如25）時，DDMRG在某些情況下誤差較大，甚至達(dá)到了50%，相比之下，解析線性響應(yīng)DMRG的精度遠(yuǎn)遠(yuǎn)大于DDMRG，隨著M的增大，兩種方法都逐漸收斂，當(dāng)M=250時，DDMRG的結(jié)果略優(yōu)于解析線性響應(yīng)DMRG.此外，Chan等［14］將DDMRG應(yīng)用到一維氫鏈模型以及Hubbard模型，通過求解單電子格林函數(shù)得到光電子能譜，以FCI的結(jié)果為基準(zhǔn)比較了DDMRG與TDDMRG的表現(xiàn)，如圖12所示，隨著M的增大，兩者均逐漸收斂到FCI的結(jié)果，DDMRG明顯收斂速度更快，此外，他們還將DDMRG用于精確計算水分子中O1s軌道的核電離能.

Fig.12 Dependence ofthe local density of states at the central site of a hydrogen chain(N=10)with different bond dimensions[14]

5.2 電聲子耦合問題

對于π電子共軛的有機分子聚集體系的能量轉(zhuǎn)移與電荷轉(zhuǎn)移常常伴隨著核的運動，其光電響應(yīng)性質(zhì)從而受到電聲子耦合的影響，當(dāng)電聲子耦合比較弱時，基態(tài)可以考慮為在被聲子云“綴飾”的接近自由的電子的圖像，當(dāng)電聲子耦合比較強時，晶格畸變導(dǎo)致電子自陷而形成了一種準(zhǔn)粒子，也被稱為極化子.White等［43］基于Lanczos DMRG研究了Holstein模型中電聲耦合強度對于電子運動的影響（圖13），圖中γ為第i個格點對應(yīng)的電聲耦合強度，t為最近鄰格點的跳躍積分.該圖展示了與光學(xué)電導(dǎo)率相關(guān)的部分物理量隨電聲耦合強度γ的變化.對于光學(xué)電導(dǎo)率有σ(ω)=Dδ(ω)+σ′(ω)，D來源于光電導(dǎo)率的Drude峰，對應(yīng)于光電導(dǎo)率相干項的權(quán)重，σ′(ω)對應(yīng)于非相干項的權(quán)重.在圖13（A）中，T為每個格點上的平均電子動能，并與σ(ω)存在的關(guān)系，即T刻畫了光電導(dǎo)率的總權(quán)重.當(dāng)電聲耦合強度γ=0時，非相干項σ′(ω)對于光電導(dǎo)率的貢獻為0，存在D=T的關(guān)系，即對應(yīng)于無相互作用的自由電子系統(tǒng).隨著γ的增加，D與T均下降，D下降得相對更快，且在γ＞2t時變得很小，光電導(dǎo)率的相干項的貢獻（D與T的比值）逐漸減小，即非相干項比重增大，由于電聲相互作用的增強形成了局域的極化子.從圖13（B）也可見，在極化子范圍內(nèi)（γ=2.5t），譜的形狀更為復(fù)雜.

Fig.13 Dynamical properties of Holstein model[43]

最近，我們［22］基于MPS/MPO的框架給出了零溫度及有限溫度下的DDMRG算法，并通過計算偶極-偶極關(guān)聯(lián)函數(shù)得到分子聚集體的吸收、發(fā)射光譜，該算法已在DMRG開源程序Renormalizer［58］中實現(xiàn).Holstein模型被用以描述同時存在分子間激子耦合與分子內(nèi)電子-振動相互作用的體系如下：

式中，εi為第i個分子的絕熱激發(fā)能；Jij是第i，j個分子之間的激子耦合；ωin與gin分別為第i個分子的第n個諧振子振動模式的頻率以及對應(yīng)的電聲耦合系數(shù)，如圖14所示.

圖15給出了DDMRG計算分子二聚體的吸收、發(fā)射強度在不同的參數(shù)空間的誤差（以精確對角化結(jié)果為標(biāo)準(zhǔn)），并與此前發(fā)展的高精度含時密度矩陣重正化群（TD-DMRG方法）［15］進行了比較.

該二聚體模型體系帶有一個頻率為ω0的簡諧振動模式（16個振動能級），選取了不同的激子耦合J=±ω0（H/J聚集體），不同的黃昆因子S=g2以及不同的溫度kBT∈[ω0，2ω0，4ω0]進行研究，對于DDMRG及TD-DMRG的計算采用了同樣的Lorentz展寬η=0.1ω0.在給定M=120的條件下，DDMRG的誤差在整體的參數(shù)空間均小于TD-DMRG，尤其是在低溫區(qū)達(dá)到了近乎精確的結(jié)果，但其誤差隨著溫度升高而上升，這是由于在高溫下存在更多的態(tài)之間躍遷，這對于CV的截斷是不利的，需要保留更大的維數(shù)（M）來達(dá)到更高的精度，因此DDMRG更適合于計算低溫區(qū)的響應(yīng)性質(zhì).需要指出的是，對于大量的有機共軛分子存在著ω0在1400 cm-1的振動序列（Vibronic progression），在室溫下處于kBT≤ω0的區(qū)間.在此之前，n粒子近似作為一種經(jīng)濟有效的方法被成功應(yīng)用于一系列光譜性質(zhì)的理論研究［59］，將DDMRG與n粒子近似進行了比較，以揭示n粒子近似的適用區(qū)間.圖16比較了雙粒子近似（2-PA）與DDMRG在3種不同的激子耦合強度下的不同溫度下的0-0發(fā)射峰的強度，可見雙粒子近似對于相對弱的耦合與較低溫度下可以給出令人滿意的結(jié)果，但在較高溫度與強激子耦合條件下表現(xiàn)較差，只考慮雙粒子近似不足以描述較大的極化子.進一步探究了0-0發(fā)射峰與激子相干長度的關(guān)系，發(fā)現(xiàn)在3種不同耦合情況二者都幾乎呈線性關(guān)系.針對可進行精確對角化的五聚體體系（最近鄰作用），比較了n粒子近似方法與DDMRG的誤差，由圖17（A）所示，當(dāng)kBT=ω0時，雙粒子近似約存在30%的高估，最簡單的單粒子近似計算在kBT=ω0時具有相當(dāng)大的誤差，使用M=50，DDMRG便實現(xiàn)了與4粒子近似相當(dāng)?shù)木?，此外，圖17（B）展示了kB=ω0時的發(fā)射強度.

Fig.14 Graphical representation of molecular aggregates characterized by Holstein model

Fig.15 Relative error for J-and H-aggregates with different Huang-Rhys factor（S∈[1.0，3.0，5.0]）across different temperature（kBT∈[0.5ω0，ω0，2ω0]）［22］

Fig.16 0-0 emission strength I0-0 as a function of T-1/2 for different excitonic coupling J0(A linear relation between them is predicted by strong excitonic coupling perturbation theory[59])[22]

作為一種新的頻率空間的響應(yīng)計算方法，CheMPS在電子-振動耦合體系的應(yīng)用仍未曾有探索，在這里，我們首次將CheMPS應(yīng)用于計算Holstein模型在零溫度以及有限溫度下的吸收光譜，展示該算法的主要特點，并做出相應(yīng)的討論.以每個分子帶有一個頻率為ω0的振動模式（振動能級數(shù)為16）的二聚體為研究對象，圖18給出了利用CheMPS在系統(tǒng)總的能量空間下計算的零溫度以及有限溫度（kBT=ω0/2）時的吸收光譜S(ω)［式（15）和式（53）］.為了將其與DDMRG進行比較，首先利用Lorentz系數(shù)［式（33），λ=4］對于給定的Chebyshev矩的個數(shù)N（即演化的步數(shù)）得到其在頻率ω處相應(yīng)的Lorentz展寬（詳見3.4節(jié)），而后進行在該展寬下的DDMRG單個頻率計算.CheMPS與DDMRG在零溫度以及有限溫度在同樣的M下實現(xiàn)了很好的吻合（且兩種方法的結(jié)果均關(guān)于M收斂）.圖19展示了不同的演化步數(shù)N得到的吸收光譜，根據(jù)在3.3節(jié)中提到的Lorentz展寬演化的步數(shù)N的關(guān)系（展寬與演化的步數(shù)呈反比），CheMPS方法可以通過調(diào)整演化的步數(shù)精確地調(diào)控光譜的精細(xì)程度.

Fig.17 Comparison between DDMRG and n-particle approximation［22］

Fig.18 Linear absorption spectrum（omit the pre-factor of the frequency index）of dimer model using CheMPS and the comparison with DDMRG

圖18所示的光譜所在的頻率區(qū)間只占了系統(tǒng)能量空間的一小部分，由圖19可見，在相當(dāng)大的能量范圍內(nèi)是沒有吸收強度的.根據(jù)3.3節(jié)提到的Lorentz展寬與實際計算時被投影的頻率區(qū)間寬度的關(guān)系，寬度越大，為了達(dá)到相同的Lorentz展寬就需要演化更多的步數(shù).對于我們在本文研究的Holstein模型而言，系統(tǒng)的能量區(qū)間的最大值是激發(fā)態(tài)的能量最高點與基態(tài)的能量最低點的差值，在模式較多或振動能級數(shù)較高的時候，如此高能態(tài)的振子強度幾乎為零，我們實際關(guān)心的光譜范圍將遠(yuǎn)遠(yuǎn)小于系統(tǒng)總的能量區(qū)間.因此，不將系統(tǒng)的整個能量區(qū)間投影到［-1，1］之間，而是選取一個有效的能量區(qū)間W*（W*＜W），尤其是對于考慮電子-振動耦合系統(tǒng)的光譜將大大減少所需演化的步數(shù).

Fig.19 Linear absorption spectrum using different numbers of Chebyshev moments for expansion

圖20展示了當(dāng)選取的有效能量區(qū)間W*為系統(tǒng)的總能量區(qū)間W的一半時的零溫度的吸收光譜，并將其與選擇總能量區(qū)間W直接計算的結(jié)果進行了比較.當(dāng)選擇的有效能量區(qū)間W*＜W時，確實可以用較少的演化步數(shù)得到同樣的精細(xì)程度（見圖20的藍(lán)線）.然而，對于同樣選取W*=W/2，如果不做能量截斷，在演化到N=70時便發(fā)散了（見圖20插圖）.隨著演化的繼續(xù)，選取有效能量空間W*的方案會產(chǎn)生誤差（見圖20的紅線），這與Krylov子空間的維度dk與掃描圈數(shù)ns有關(guān)，圖21展示了不同的dk與ns計算光譜的相對誤差，可以發(fā)現(xiàn)，dk與ns對于誤差的影響均不是單調(diào)的，因為截斷是局部進行的，因此可能會存在“欠截斷”與“過截斷”的效果.說明選取怎樣的有效能量空間W*，多大的dk與ns是沒有明確的準(zhǔn)則去遵循的，而且其最優(yōu)組合將視不同體系決定.

Fig.20 Linear absorption spectrum calculated by using effective frequency window and the energy truncation scheme with the Krylov subspace dimension dk=10 and ns=1 sweep times

Fig.21 Relative error using effective band(W*=W/2)at N=700 with respect to the full many body band calculationat N=1000 for different combinations of Krylov subspace dimension dk and sweep times ns for energy truncation

6 總結(jié)與展望

DMRG作為一種計算低維強關(guān)聯(lián)體系電子結(jié)構(gòu)的方法，近20多年來已經(jīng)逐漸地成為量子化學(xué)的重要計算手段.對于電子結(jié)構(gòu)問題，DMRG可以作為CASSCF中取代Full CI的求解器，從而使得可處理的活性空間遠(yuǎn)大于傳統(tǒng)的方法.近幾年，隨著含時DMRG和響應(yīng)理論的發(fā)展，已迅速發(fā)展成為計算復(fù)雜體系量子動力學(xué)和光譜的高精度方法.求解響應(yīng)性質(zhì)有兩種不同的角度：（1）直接在頻率域求解響應(yīng)函數(shù)；（2）通過時間演化得到時間關(guān)聯(lián)函數(shù)，并通過傅里葉變換得到頻率域的信息.本文詳細(xì)闡述了頻率空間的幾種不同的算法，并介紹了其應(yīng)用.DDMRG的優(yōu)點是可以實現(xiàn)大規(guī)模的并行計算，即每個頻率的精確計算都是獨立進行的，最后收集起來就得到了整個頻域的譜線.對于只需計算少量點的響應(yīng)性質(zhì)計算的時候，DDMRG將是首要的選擇.CheMPS盡管也是一種頻率域的計算方法，但其可以一次性得到所有頻率的響應(yīng)信息，因其更為快速，將CheMPS首次應(yīng)用到計算電聲子耦合系統(tǒng)的響應(yīng)性質(zhì)，展示了該算法的主要特點以及進一步探索的方向，如果系統(tǒng)有響應(yīng)信息的頻率區(qū)間占整個系統(tǒng)的能量空間的絕大部分，CheMPS將是一種兼顧精度與計算量的計算方法，如果有響應(yīng)的頻率區(qū)間只占了整個系統(tǒng)的能量空間的一小部分，可以在有效的頻率空間進行展開，同時需要做能量截斷，這可以大大節(jié)省演化的步數(shù)，然而一種高效的、可靠的能量截斷方案有待進一步的探索.電聲耦合體系擁有較大的系統(tǒng)能量空間，尤其是對于多分子、多模式與振動能級高的情況，需要探索一種高效的、可靠的能量截斷方案.

TD-DMRG作為一種數(shù)值精確的方法被廣泛用于超快動力學(xué)的計算，相較于頻率空間的方法，TD-DMRG可以給出實時動力學(xué)信息.此外，TD-DMRG可以一次性得到整個頻率空間的響應(yīng)性質(zhì).TD-DMRG的弱點在于長時間演化，首先在演化的過程中存在誤差的累積，另外，由于系統(tǒng)的糾纏隨著時間增加，為了獲得相同的精度，所需的M應(yīng)隨時間增加而變大［33～35］.頻率空間的方法可以避免這一問題，有望給出更為精確的結(jié)果.在頻率空間算法中，CV-DMRG以及其衍生算法DDMRG算法需要對每一個頻率分別進行計算，具有非常高的精度，但是存在計算速度慢的問題，但可以利用其天然并行的優(yōu)勢以及通過GPU加速張量計算［22］大大縮短計算時間.

借助于熱場動力學(xué)（純化）的方法，響應(yīng)性質(zhì)的計算被自然地推廣到有限溫度.熱場動力學(xué)方法以無窮高溫度下的最大糾纏態(tài)作為演化的初態(tài)進行虛時演化，因此對于較低溫度需要演化更多步，最小糾纏典型量子熱態(tài)（Minimally entangled typical thermal states，METTS）方法［60］是對于熱場動力學(xué)方法在T→0的低溫計算的有力補充.

近幾年的發(fā)展也表明，DMRG方法尤其是其矩陣乘積態(tài)的理論形式除了作為一種解決電子相關(guān)問題的有效方法，也是一種高效、準(zhǔn)確的復(fù)雜體系量子動力學(xué)和光譜計算方法，這方面的研究還處于初始階段，大量的探索性工作有待開展，本文所涉及頻域的動態(tài)DMRG代表著一種探索，將在更多的應(yīng)用實例中得到發(fā)展.