李維唐，任佳駿，帥志剛

（清華大學(xué)化學(xué)系,有機(jī)光電子與分子工程教育部重點(diǎn)實(shí)驗(yàn)室,北京100084）

1992年，由White［1，2］提出的密度矩陣重正化群（Density matrix renormalization group，DMRG）算法是研究一維量子多體強(qiáng)關(guān)聯(lián)系統(tǒng)的強(qiáng)有力方法，尤其適用于只有最近鄰相互作用的量子模型.隨后，Shuai，Ramasesha和Fano等將DMRG推廣到具有長(zhǎng)程作用勢(shì)的化學(xué)體系［3～6］.通過(guò)將DMRG的基矢對(duì)稱化，Shuai等［3～5］計(jì)算了共軛多烯的低激發(fā)態(tài)結(jié)構(gòu)，為高分子發(fā)光給出了判據(jù)，并發(fā)展了基于DMRG的線性和非線性光學(xué)響應(yīng)系數(shù)的計(jì)算方法.1999年，White和Martin［7］首次將DMRG用于從頭算量子化學(xué)模型，隨后Head-Gordon，Chan，Reiher，Yanai和Kurashige等［8～11］對(duì)DMRG算法做了大量的優(yōu)化，使得量子化學(xué)DMRG成為了強(qiáng)關(guān)聯(lián)分子體系的一種主流的方法.含時(shí)DMRG（Time-dependent DMRG，TD-DMRG）是針對(duì)分子器件的量子輸運(yùn)現(xiàn)象而提出［12］，并得到快速發(fā)展［13，14］.早期的TD-DMRG時(shí)間演化算法如時(shí)間步瞄準(zhǔn)（Time-step targeting，TST）算法主要是受最初的DMRG迭代算法啟發(fā)［15］，后期發(fā)展的大多數(shù)TD-DMRG算法是基于矩陣乘積態(tài)（Matrix product state，MPS）的波函數(shù)擬設(shè)和矩陣乘積算符（Matrix product operator，MPO）設(shè)計(jì)的［16，17］，包括一系列演化再壓縮（Propagation and compression，P&C）算法［18］以及最近的基于含時(shí)變分原理（Time-dependent variational principle，TDVP）算法［19，20］.只要MPS的虛擬鍵維數(shù)足夠大，這些TD-DMRG算法可以幾乎精確地給出復(fù)雜多體系統(tǒng)的含時(shí)薛定諤方程的解.

近年來(lái)，TD-DMRG主要被用于研究電子-聲子（振動(dòng)）耦合問(wèn)題（簡(jiǎn)稱電聲耦合）.電聲耦合在化學(xué)反應(yīng)、（生物）分子的光化學(xué)光物理過(guò)程以及固體中的輸運(yùn)現(xiàn)象中均扮演重要角色.如TD-DMRG已經(jīng)被成功應(yīng)用于模擬給體/受體異質(zhì)結(jié)的激子分離［21，22］，分子聚集體的熒光和吸收光譜［23，24］以及有機(jī)半導(dǎo)體的電荷遷移率［25］.實(shí)際上，由于電聲耦合問(wèn)題在化學(xué)中的地位至關(guān)重要，早在DMRG之前，理論化學(xué)家已經(jīng)開發(fā)了一整套適用于不同科學(xué)問(wèn)題的理論方法來(lái)模擬電聲耦合動(dòng)力學(xué)過(guò)程，即非絕熱動(dòng)力學(xué)方法［26，27］.這類方法的目標(biāo)是超越玻恩-奧本海默近似，對(duì)電子和原子核的耦合運(yùn)動(dòng)進(jìn)行描述.眾多的非絕熱動(dòng)力學(xué)方法大體上可以分為兩類：（1）利用量子力學(xué)處理原子核的運(yùn)動(dòng)，稱為量子動(dòng)力學(xué)方法，其優(yōu)點(diǎn)是可以正確描述隧穿效應(yīng)及零點(diǎn)振動(dòng)，代表性方法是多組態(tài)含時(shí)Hartree（Multi-configuration time-dependent Hartree，MCTDH）［28，29］及其多層推廣（Multi-layer MCTDH，ML-MCTDH）［30，31］；（2）利用經(jīng)典力學(xué)近似描述原子核的運(yùn)動(dòng)，稱為混合量子/經(jīng)典動(dòng)力學(xué)方法，其優(yōu)點(diǎn)是計(jì)算量較小，代表性方法是勢(shì)能面跳躍法［32，33］.TD-DMRG利用波函數(shù)描述原子核的狀態(tài)，因此屬于量子動(dòng)力學(xué)方法.TD-DMRG與（ML-）MCTDH在非絕熱動(dòng)力學(xué)中的定位類似，其理論形式也有共通之處［34］.在最近幾年，TD-DMRG與（ML-）MCTDH相互促進(jìn)，產(chǎn)生了一批新的算法［35～37］.

本文將簡(jiǎn)述基于矩陣乘積態(tài)和矩陣乘積算符的DMRG的基礎(chǔ)理論，并將分3類較詳細(xì)地介紹現(xiàn)有的TD-DMRG時(shí)間演化算法，然后，介紹在TD-DMRG中引入溫度效應(yīng)的算法，最后，討論TD-DMRG算法的應(yīng)用，包括我們［23，25］最近利用TD-DMRG模擬分子聚集體光譜和有機(jī)半導(dǎo)體遷移率的一些進(jìn)展.本文聚焦于TD-DMRG在電聲耦合問(wèn)題中的應(yīng)用，如關(guān)注溫度效應(yīng)是因?yàn)闇囟葘?duì)振動(dòng)自由度有重要影響.然而需要強(qiáng)調(diào)的是，不能簡(jiǎn)單地認(rèn)為TD-DMRG就是一種新興的可以作為（ML-）MCTDH替代品的非絕熱動(dòng)力學(xué)方法，因?yàn)門D-DMRG算法本身并不限制于特定的模型.如最近化學(xué)領(lǐng)域的一些研究顯示出TD-DMRG模擬基于第一性原理哈密頓的電子動(dòng)力學(xué)的潛力［38，39］.盡管這些研究目前只占使用TD-DMRG在化學(xué)中應(yīng)用的少數(shù)，但其應(yīng)用領(lǐng)域?qū)?huì)持續(xù)拓寬.關(guān)于TD-DMRG的理論與應(yīng)用的更詳細(xì)內(nèi)容可參見相關(guān)報(bào)道［17，40］.

1 基于MPS的現(xiàn)代DMRG理論

現(xiàn)代DMRG理論一般基于矩陣乘積態(tài)擬設(shè)［16］.對(duì)于由N個(gè)自由度組成，每個(gè)自由度的正交歸一基為的系統(tǒng)，任意量子態(tài)Ψ可以表示成N個(gè)矩陣的連乘積，即矩陣乘積態(tài)：

式中：{·}表示各個(gè)指標(biāo)的集合.圖1為矩陣乘積態(tài)的圖示.ai的維數(shù)被稱為（虛擬）鍵維數(shù)，記為M或|ai|.

Fig.1 Graphical representation of an MPS

多體重正化基定義為

對(duì)于一個(gè)特定的量子態(tài)而言，其矩陣乘積態(tài)表示并不是唯一的.在任意兩個(gè)相鄰矩陣AσiAσi+1之間插入一個(gè)單位矩陣I=GG-1不改變波函數(shù)，但局域矩陣卻發(fā)生了變化：所張成的空間的系數(shù)矩陣.通常，

為方便計(jì)算，可以引入額外的規(guī)范條件來(lái)消除參數(shù)的冗余性.最常見的規(guī)范條件是“混合/左/右規(guī)整化”.一個(gè)規(guī)整中心在第n個(gè)矩陣的混合規(guī)整矩陣乘積態(tài)可以寫為

Lσi(Rσj)*表示Lσi(Rσj)的復(fù)共軛.

當(dāng)滿足式（8）和式（9）時(shí)，重正化基滿足正交歸一條件S[1：i]lil′i=δlil′i(i=1，2，…，n-1）和S[j+1：N]rjr′j=δrjr′j(j=n，n+1，…，N-1）.當(dāng)規(guī)整中心n處于最右（左）側(cè)時(shí)，該MPS被稱為是左（右）規(guī)整的.若Cσn進(jìn)一步被QR分解：

式中：Lσn滿足式（8），而Dlnrn是空間的系數(shù)矩陣.之后，Dlnrn與相乘，得到，即規(guī)整中心的位置向右移動(dòng)了一個(gè)位置.通過(guò)與上述過(guò)程相似的RQ分解，可以將規(guī)整中心的位置向左移動(dòng)一個(gè)矩陣（以下將QR分解和RQ分解統(tǒng)稱為QR分解）.通常，一個(gè)規(guī)整中心在第n個(gè)矩陣上的MPS可以由任意MPS做連續(xù)的QR分解得到，這一過(guò)程稱為規(guī)整化.

若利用奇異值分解（Singular value decomposition，SVD）代替式（10）中的QR分解：

式中：Λ是一個(gè)對(duì)角元從大到小排列的實(shí)對(duì)角矩陣并且若是歸一化的，則滿足式（8）.若所有的都被保留，并且向右移動(dòng)規(guī)整中心那么MPS所表示的波函數(shù)不發(fā)生變化，SVD的作用與QR分解相同.然而，如果只保留最大的M項(xiàng)（M＜k）以得到那么從數(shù)學(xué)上來(lái)說(shuō)將會(huì)是的一個(gè)好的近似，且M更小.若是歸一化的，由截?cái)嗨鶐?lái)的誤差可以由衡量.基于這個(gè)算法，一個(gè)左（右）規(guī)整的矩陣乘積態(tài)可以通過(guò)從第N個(gè)到第1個(gè)（從第1個(gè)到第N個(gè)）矩陣的連續(xù)SVD“壓縮”得到其M更小的近似在每一步截?cái)鄷r(shí)，有兩種常用的標(biāo)準(zhǔn).第一種是截?cái)嘀令A(yù)先指定的虛擬鍵維數(shù)M.第二種是保留所有大于預(yù)先指定閾值的Λss.

除了規(guī)整化操作和壓縮操作，MPS結(jié)構(gòu)還允許其它一些常見操作，如與MPS相似，常見的量子算符可以表示成矩陣乘積算符（MPO）的形式［16，41］：

圖2為MPO的圖示.

Fig.2 Graphical representation of an MPO

其中

其中

2 時(shí)間演化算法

TD-DMRG的時(shí)間演化算法種類繁多，但大體上可以分為3類.（1）對(duì)時(shí)間演化算符e-iHt或進(jìn)行基于MPO/MPS全局的近似.這類方法包括龍格-庫(kù)塔（Runge-Kutta，RK）方法［18，23］、時(shí)間演化塊消減（Time-evolving block decimation，TEBD）方法［14，42，43］、WI，II方法［44］、Krylov子空間方法［18，45］、Chebyshev展開方法［46］以及分裂算符方法等［47］.這些方法的共通點(diǎn)是在每一步時(shí)間演化時(shí)，波函數(shù)MPS首先通過(guò)MPO/MPS乘法和MPS/MPS加法作為一個(gè)整體演化（Propagate）形成具有更大虛擬鍵維數(shù)的MPS，隨后依據(jù)原本的M或奇異值截?cái)嚅撝颠M(jìn)行壓縮（Compress）.因此這類算法被稱為P&C算法.（2）基于含時(shí)變分原理（TDVP）［48］.依據(jù)不同的推導(dǎo)運(yùn)動(dòng)方程的方式，這類方法有固定規(guī)范自由度［19］和投影算符裂分（Projector-splitting，PS）［20］兩個(gè)變種.（3）受原本的基態(tài)DMRG算法啟發(fā)，利用平均約化密度矩陣來(lái)更新重正化基.其代表性方法是時(shí)間步瞄準(zhǔn)方法［15］及一些變種［38，49］.除了TEBD之外，所有上述時(shí)間演化算法都適用于具有長(zhǎng)程相互作用的模型.此外P&C算法是在現(xiàn)代MPS/MPO框架下最直接的時(shí)間演化算法，而PS算法是最近最常被使用的算法［24，25，50～52］.從數(shù)值計(jì)算的角度上來(lái)說(shuō)，所有的時(shí)間演化算法都需要進(jìn)行大量的矩陣乘法運(yùn)算，可以有效被圖形處理器（Graphical processing units，GPU）所加速.GPU加速不同類型算法的加速比有差異.不需要壓縮維數(shù)較大的MPS的算法可以最大程度上被GPU加速，加速比可以達(dá)到數(shù)十倍［37］.這類算法主要是基于TDVP的算法.

2.1 演化并壓縮算法

為簡(jiǎn)便起見，采用原子單位，設(shè)?為1.基于最簡(jiǎn)單的前向歐拉法，對(duì)每一時(shí)間步τ應(yīng)該按照下式演化Ψ：

正如前文所介紹，MPO/MPS乘法和MPS/MPS加法是MPS的常規(guī)操作，因此式（18）可以很方便地轉(zhuǎn)換成一個(gè)TD-DMRG算法.通過(guò)式（18）獲得的往往比初始態(tài)具有更大的虛擬鍵維數(shù)，所以在下一步迭代之前需要對(duì)進(jìn)行壓縮.這一具有兩個(gè)步驟的算法因而被稱為是P&C算法.

對(duì)這一簡(jiǎn)單算法的一個(gè)直接改進(jìn)是將具有一階誤差的前向歐拉積分法替換為經(jīng)典的具有四階誤差的龍格-庫(kù)塔（RK4）積分法［18，23］：

對(duì)于不含時(shí)的哈密頓來(lái)說(shuō)，RK4算法與形式時(shí)間演化算符e-iHτ的四階泰勒展開等價(jià)：

在實(shí)際計(jì)算時(shí)，第n階項(xiàng)可以由第n-1階項(xiàng)作用H計(jì)算得到.

其它P&C算法雖然理論基礎(chǔ)各不相同，但具體到TD-DMRG算法與本文介紹的基本類似.其中，TEBD算法是最著名的P&C算法之一［14，42，43］，雖然該算法具有很高的效率且易于實(shí)現(xiàn)，但是其不能被直接用于具有長(zhǎng)程相互作用的模型，因此在化學(xué)中的應(yīng)用受到一定限制［21，22］.假設(shè)一個(gè)具有2N個(gè)自由度的系統(tǒng)可以由只有最近鄰相互作用的哈密頓描述：

式中：hj，j+1作用于第j和第j+1個(gè)自由度.觀察到哈密頓可以分為兩個(gè)部分H=H1+H2：

按照這種劃分方式，一階Trotter分解可以將e-iHτ分為兩部分乘積：

盡管H1和H2相互并不對(duì)易，H1和H2中的各個(gè)項(xiàng)相互對(duì)易，因此：

連乘中的每一項(xiàng)e-ihτ都可以表示成一個(gè)兩個(gè)格點(diǎn)的MPO，即都可以被有效地表示為MPO.將這些算符作用到上后進(jìn)行壓縮就完成了時(shí)間演化中的一步.同理也可以構(gòu)造二階Trotter分解：

對(duì)于具有長(zhǎng)程相互作用的體系，TEBD需要和其它技術(shù)結(jié)合，如引入交換門［53］或在電聲耦合模型中對(duì)系統(tǒng)哈密頓做變換［54］.此外，Trotter分解的技術(shù)也可以和泰勒展開時(shí)間演化算符的做法相結(jié)合［47］.其思路是將動(dòng)能部分從哈密頓中分離出去，以減小時(shí)間演化算符中的MPO的虛擬鍵維數(shù).

WI，II方法［44］試圖基于哈密頓H的形式顯式地將時(shí)間演化算符e-iHτ近似為一個(gè)MPO.該方法與RK算法相比的優(yōu)勢(shì)是其每一個(gè)自由度平均的積分誤差是一個(gè)常數(shù).

Krylov子空間方法（也被稱為L(zhǎng)anczos方法）是近似的有效方法［18，45］.Krylov子空間的定義是由向量所張成的空間，維數(shù)是N.這些向量一般不是正交歸一的，正交歸一化后得到的向量稱為Krylov向量{k0，k1，…，kN-1}.Krylov子空間方法的總體目標(biāo)是找到在Krylov子空間中的最佳近似.當(dāng)N等于原本的希爾伯特空間維數(shù)時(shí)，該方法就是精確的.在實(shí)際計(jì)算中，往往少數(shù)Krylov向量就可以得到誤差非常小的結(jié)果.向Krylov子空間投影的投影算符定義如下：

Chebyshev展開方法利用Chebyshev多項(xiàng)式展開e-iHt［46］.Chebyshev向量依據(jù)以下遞推關(guān)系計(jì)算：

式中，φ0(t)=c0(t)，φn＞0=2cn(t)，而cn(t)的定義為

式中，Jn是第一類Bessel函數(shù).

在所有P&C時(shí)間演化方法中，目前應(yīng)用最廣泛的是基于e-iHτ的泰勒展開的方法，TEBD方法和Krylov子空間方法.對(duì)于只具有最近鄰相互作用或可以轉(zhuǎn)化成只具有最近鄰相互作用的哈密頓來(lái)說(shuō)，TEBD在易于實(shí)現(xiàn)和較高精度間取得了較好的平衡［51，55～57］.對(duì)于一般的哈密頓而言，基于e-iHτ的泰勒展開的方法最容易實(shí)現(xiàn)［23，47］，但精度不如Krylov子空間方法［39］.

2.2 基于含時(shí)變分原理的算法

Rayleigh-Ritz變分原理被廣泛用于尋找不含時(shí)哈密頓的近似基態(tài).類似地，含時(shí)變分原理（TDVP）也是在已知初始態(tài)時(shí)尋找最佳的近似含時(shí)波函數(shù)的有力工具.Dirac-Frenkel含時(shí)變分原理［48，58］為

在實(shí)時(shí)演化中，TDVP可以保證波函數(shù)的模和體系能量守恒［58］.這兩個(gè)特性對(duì)于可靠的長(zhǎng)時(shí)間動(dòng)力學(xué)非常重要.從幾何的角度而言，TDVP可以理解成是將正交投影到在當(dāng)前時(shí)間的切空間中：

其中，P是由切空間中正交歸一的向量構(gòu)造成的投影算符，定義為

圖3為式（33）的圖形表示.重疊矩陣的逆S-1是由重正化基的非正交性導(dǎo)致的，而減號(hào)“-”里的項(xiàng)是為了消除參數(shù)冗余性［19，59］.

Fig.3 Graphical representation of the tangent space projector

兩種基于TDVP的時(shí)間演化方式區(qū)別在于在解式（32）時(shí)采用了不同的MPS規(guī)范條件.在第1種時(shí)間演化方案中，MPS的規(guī)范條件不發(fā)生改變.為方便起見，可以將除i=n的一個(gè)“+”項(xiàng)之外其它相鄰“+”項(xiàng)和“-”項(xiàng)合并，如此每一項(xiàng)都是正交的，將式（33）中的投影算符轉(zhuǎn)變?yōu)?/p>

其中

這種形式的投影算符及其相應(yīng)的切空間向量首先由Haegeman等［19］針對(duì)均一MPS（Uniform MPS，uMPS）提出，Wouters等［59］對(duì)其進(jìn)行了重新描述以推導(dǎo)一些針對(duì)基態(tài)和激發(fā)態(tài)的Post-DMRG算法.

如果采取一定的規(guī)范條件，式（39）和式（40）中的一些重疊矩陣將變成單位陣，表達(dá)式可以得到簡(jiǎn)化.假設(shè)當(dāng)前MPS是左規(guī)整的，規(guī)整中心在第N個(gè)矩陣處，則S[1：i]簡(jiǎn)化為I，且在式（38）中設(shè)n=N最為簡(jiǎn)便.將簡(jiǎn)化后的投影算符代入式（32）得：

其中

式（41）和式（42）構(gòu)成了一組耦合的非線性方程，與（ML-）MCTDH中標(biāo)準(zhǔn)的運(yùn)動(dòng)方程十分類似［28～30］.（ML-）MCTDH中有兩種典型的時(shí)間演化方案，一種稱為可變平均場(chǎng)（Variable mean field，VMF），另一種稱為常數(shù)平均場(chǎng)（Constant mean field，CMF）.這兩種方案也可用于對(duì)式（41）和式（42）進(jìn)行積分［37］.VMF使用普適的初值算法如前向歐拉法直接求解這些耦合方程.TDVP首先應(yīng)用于MPS時(shí)即采用的這種方案［19］.CMF則假設(shè)H[i]和S[i+1：N]比起CσN，Lσi的變化速率慢很多.因此在時(shí)間步τ中，“平均場(chǎng)”H[i]，S[i+1：N]可以認(rèn)為是常數(shù)，而演化局域矩陣CσN，Lσi.CMF因而可以被認(rèn)為是VMF的一種近似.

這種固定規(guī)范條件的演化算法中，對(duì)S求逆的運(yùn)算需要特別注意，因?yàn)槿绻鸖的一些特征值很小，對(duì)S求逆將會(huì)是數(shù)值不穩(wěn)定的.當(dāng)體系僅僅處于弱關(guān)聯(lián)狀態(tài)或虛擬鍵維數(shù)非常大時(shí)，比如體系處于Hartree直積態(tài)時(shí)，該問(wèn)題會(huì)變得比較嚴(yán)重.在某種程度上，這種數(shù)值不穩(wěn)定問(wèn)題使得這種時(shí)間演化算法處于一種悖論中：使用更大的虛擬鍵維數(shù)原則上應(yīng)該使計(jì)算結(jié)果更精確，然而如果時(shí)間步長(zhǎng)不改變的話，使用更大的虛擬鍵維數(shù)實(shí)際上會(huì)使結(jié)果變得更不精確.相同的問(wèn)題在（ML）-MCTDH中也會(huì)出現(xiàn)，為了使運(yùn)動(dòng)方程積分更穩(wěn)定，S一般由一個(gè)正則化的重疊矩陣代替［29］：

式中：ε是一個(gè)小標(biāo)量，一般取值為10-8～10-14.最近，Meyer和Wang［60，61］提出了一種基于對(duì)系數(shù)矩陣進(jìn)行SVD矩陣展開（Matrix unfolding，MU）的改進(jìn)（ML-）MCTDH正則化方式.這種正則化方式可以使時(shí)間積分算法更加精確和穩(wěn)定（Robust）.相同的方式可被用于積分式（42），這種時(shí)間演化算法稱為TDVPMU［37］.當(dāng)計(jì)算重疊矩陣時(shí)，可以將當(dāng)前MPS的一份拷貝的規(guī)范中心移動(dòng)到第i+1個(gè)矩陣，且該矩陣可進(jìn)一步被SVD分解：

其中，“[…]”中的表達(dá)式正是S[i+1：N]-1.這種新的正則化方式的關(guān)鍵點(diǎn)在于具有下劃線的部分可以

此處的1/2次方是為了與式（48）中原本的正則化方式保持一致.式（50）中的則保持不變，以減少對(duì)運(yùn)動(dòng)方程的影響，與MCTDH文獻(xiàn)中的做法一致［60，61］.盡管在MU中需要對(duì)環(huán)境做規(guī)整化，進(jìn)行時(shí)間演化的MPS的規(guī)范條件保持不變.

第2種基于TDVP的時(shí)間演化方案是投影算符裂分法（PS）.投影算符裂分的基礎(chǔ)在于式（33）中的各項(xiàng)投影算符可以采取不同的規(guī)范條件.在對(duì)一個(gè)一般的MPS做從第N個(gè)矩陣到第i+1個(gè)矩陣的規(guī)整化后成為了右側(cè)的重正化的基，與的關(guān)系如下：

通過(guò)QR分解得到的矩陣D是上三角矩陣.因此，重疊矩陣S[i+1：N]等于D*DT，而式（3）中對(duì)于一個(gè)一般的MPS定義的投影算符P[i+1：N]變形為

對(duì)于P[1：i]同理可得：

如此定義的切空間投影算符不含有求重疊矩陣逆的運(yùn)算，似是對(duì)式（34）和式（35）的重大改進(jìn).然而，因?yàn)樵谶@種定義方式中規(guī)范條件對(duì)于不同項(xiàng)的投影算符是不同的，上述介紹的VMF和CMF積分方式不再適用.Lubich和Haegeman等建議使用對(duì)稱二階Trotter分解將時(shí)間演化算符裂分成單獨(dú)的項(xiàng)［20，62］：

基于式（55）中的時(shí)間演化算符，一步時(shí)間演化包括了一次從左向右的掃描以及隨后的從右向左的掃描，兩次掃描的時(shí)間步長(zhǎng)均是τ/2.以從左向右的掃描為例，位于規(guī)整中心的矩陣首先通過(guò)作用進(jìn)行正向時(shí)間演化：

其中，H[i]及其組成部分具有與式（43），式（44），式（45），式（46）相同的定義，只是式（46）中的Aσi被替換成為L(zhǎng)σi或者Rσi.然后，時(shí)間演化后的矩陣通過(guò)QR分解得到左規(guī)整的矩陣和系數(shù)矩陣通過(guò)作用投影算符P[1：i]?P[i+1：N]進(jìn)行逆向時(shí)間演化：

與搜索基態(tài)的雙格點(diǎn)DMRG算法類似，TDVP-PS也可以采用雙格點(diǎn)的算法，使得虛擬鍵維數(shù)可以隨體系糾纏熵的增長(zhǎng)而增長(zhǎng).雙格點(diǎn)算法在數(shù)值上也比單格點(diǎn)的算法更加穩(wěn)定［17，51］.但是，與基態(tài)算法類似，雙格點(diǎn)算法比單格點(diǎn)算法計(jì)算量大很多.

TDVP-MU和TDVP-PS都是基于含時(shí)變分原理，當(dāng)不考慮數(shù)值誤差時(shí)，兩者給出的時(shí)間演化結(jié)果應(yīng)該是一樣的.與P&C類算法相比，TDVP-MU以及單格點(diǎn)的TDVP-PS都需要提前指定一個(gè)固定的虛擬鍵維數(shù)，且如果當(dāng)初態(tài)關(guān)聯(lián)很弱時(shí)，需要構(gòu)建額外的重正化基加入MPS.雙格點(diǎn)的TDVP-PS一般不需要構(gòu)建額外的重正化基使虛擬鍵維數(shù)達(dá)到指定的數(shù)值，因?yàn)樵谶M(jìn)行每?jī)蓚€(gè)格點(diǎn)的時(shí)間演化時(shí)需要對(duì)代表兩個(gè)格點(diǎn)的矩陣做SVD，這就允許動(dòng)態(tài)調(diào)整這兩個(gè)格點(diǎn)之間虛擬鍵維數(shù).TDVP-MU和TDVP-PS的主要區(qū)別是TDVP-MU需要進(jìn)行額外的正則化，而TDVP-PS則不需要.

2.3 時(shí)間步瞄準(zhǔn)算法

時(shí)間步瞄準(zhǔn)（TST）算法主要是受同時(shí)求解基態(tài)和少數(shù)幾個(gè)低能激發(fā)態(tài)的態(tài)平均DMRG算法啟發(fā)，其數(shù)學(xué)基礎(chǔ)不像前文介紹的算法那樣嚴(yán)格［15，17］.本文僅簡(jiǎn)要介紹其基本思想.

此時(shí)，整體波函數(shù)仍然表示體系在t時(shí)刻的狀態(tài)，僅第i個(gè)矩陣所表示的重正化基依據(jù)從t到t+τ的時(shí)間步中的狀態(tài)進(jìn)行了更新.這一步更新不是嚴(yán)格的，會(huì)引入一定的誤差.這一掃描過(guò)程不斷進(jìn)行直到收斂，而此時(shí)MPS的基被認(rèn)為是能同時(shí)很好地描述t時(shí)刻與t+τ時(shí)刻的狀態(tài).此后，在規(guī)整中心進(jìn)行真正的時(shí)間演化，完成MPS整體的一步時(shí)間演化.

3 有限溫度算法

化學(xué)問(wèn)題多與溫度相關(guān)，如在以電聲耦合模型來(lái)描述的輸運(yùn)現(xiàn)象中，溫度起到了非常重要的作用，需要從量子統(tǒng)計(jì)的混合態(tài)的角度加以考慮.在波函數(shù)框架的TD-DMRG形式中有兩種方式可以考慮溫度帶來(lái)的混合態(tài)問(wèn)題，第1種為純化（Purification）或輔助物（Ancilla）或熱場(chǎng)動(dòng)力學(xué)（Thermal field dynamics），適合高溫和常溫下的情形［63，64］.第2種為最小糾纏典型量子熱態(tài)（Minimally entangled typical thermal state，METTS）方法，主要適合低溫下的情形［53，65］.這兩種方法也可以結(jié)合使用［66］.需要指出的是，這兩種方法都不是專門為TD-DMRG或MPS設(shè)計(jì)的，而是適用于各種基于波函數(shù)的方法.對(duì)于化學(xué)問(wèn)題，純化方法是最常用的一種直接了當(dāng)?shù)挠?jì)算方法.

3.1 在擴(kuò)增的希爾伯特空間純化混合態(tài)的算法

純化方法的核心是將混合態(tài)（密度矩陣）在擴(kuò)增的希爾伯特空間表示為純態(tài)（波函數(shù)）［63，64，67］.對(duì)于任意物理量O，其有限溫度期望的表達(dá)式為

而這正是需要用密度矩陣描述溫度效應(yīng)的理由.由于MPS是一種波函數(shù)的擬設(shè)，因此最好避免使用密度矩陣，用一個(gè)符合式（60）的波函數(shù)來(lái)描述體系性質(zhì).這一目標(biāo)可以由將輔助Q空間引入P空間得到.Q空間的結(jié)構(gòu)與P空間完全相同，可以定義相同的狀態(tài)和算符.如果Q空間的能量本征態(tài)由上波浪線來(lái)表示那么構(gòu)成了擴(kuò)增的希爾伯特空間的一組正交歸一完備基.定義如下的擴(kuò)增的希爾伯特空間中的波函數(shù)

由于O僅作用于P空間，容易證明滿足式（60）.如果需要，熱平衡態(tài)的密度矩陣可以由ΨPQ對(duì)Q求偏跡（或偏內(nèi)積）構(gòu)造：可以從P和Q空間的“單位狀態(tài)出發(fā)進(jìn)行虛時(shí)演化計(jì)算得到：在任意基組下的形式都相同，所以可采取易于計(jì)算的基組，事先并不需要計(jì)算配分函數(shù)也不需要顯式地計(jì)算出來(lái)，因?yàn)橐脖环Q為最大糾纏態(tài)，因?yàn)镻空間和Q空間之間的糾纏達(dá)到最大.由于其可以當(dāng)做是的歸一化系數(shù).在每一步虛時(shí)間演化后進(jìn)行歸一化也可以解決波函數(shù)的模在虛時(shí)演化中不守恒的問(wèn)題.在獲得后，可以對(duì)其進(jìn)行實(shí)時(shí)間演化，以獲得所需的實(shí)時(shí)性質(zhì).原則上，開展實(shí)時(shí)演化和虛時(shí)演化的順序是可以交換的，但實(shí)際上，先進(jìn)行虛時(shí)演化再實(shí)時(shí)演化的計(jì)算效率更高，此外先進(jìn)行虛時(shí)演化得到有利于對(duì)體系的熱平衡態(tài)的性質(zhì)進(jìn)行分析.

虛時(shí)演化初始態(tài)I具有最大糾纏的事實(shí)對(duì)TD-DMRG計(jì)算具有不利影響，因?yàn)榧m纏越大，MPS所需的虛擬鍵維數(shù)和計(jì)算量就越大.實(shí)際上，之所以有限溫度計(jì)算的計(jì)算量一般比零溫度計(jì)算大得多（特別是在有振動(dòng)自由度的時(shí)候），不僅是因?yàn)榧兓笞杂啥鹊臄?shù)目翻倍，更是因?yàn)樘摃r(shí)演化和實(shí)時(shí)演化所需要的虛擬鍵維數(shù)一般比零溫度大很多.幸運(yùn)的是，在化學(xué)中常出現(xiàn)的一類電聲耦合問(wèn)題中，振動(dòng)由獨(dú)立的簡(jiǎn)諧振子描述，且體系初態(tài)處于電子和振動(dòng)相互不關(guān)聯(lián)的態(tài)，其中振動(dòng)處于熱平衡態(tài).此時(shí)可以對(duì)振動(dòng)自由度的PQ空間進(jìn)行基變換，消除PQ空間之間的糾纏［68，69］.設(shè)電子自由度處于基態(tài)而體系的振動(dòng)處于熱平衡態(tài)，則體系總的密度矩陣為

式（61）定義的純化后的熱平衡態(tài)與對(duì)擴(kuò)展空間的基態(tài)進(jìn)行一個(gè)酉變換的結(jié)果相同［67］：

其中G的形式為

以化學(xué)中常見的Holstein-Peierls哈密頓為例：

經(jīng)過(guò)基變換后的哈密頓形式為

3.2 最小糾纏典型量子熱態(tài)算法

低溫下，純化方法的效率較低.在β→∞極限，一方面，熱平衡態(tài)實(shí)質(zhì)上就是基態(tài)，因此純化引入的Q空間成為了MPS表示的累贅，另一方面，β→∞本身就意味著需要做無(wú)窮多步的虛時(shí)演化［16］.最小糾纏典型量子熱態(tài)（METTS）方法是一種基于采樣的方法，可以避免純化，在低溫下效率較高［53，65］.

METTS方法的出發(fā)點(diǎn)是典型量子熱態(tài)（TTS）的概念.式（59）是有限溫度的物理量期望O在能量基下的表示.而對(duì)于任意完備正交歸一的基期望的表達(dá)式為

其中，n表示不同的自由度.這么做的原因是（盡管不嚴(yán)格）：假如虛時(shí)間演化初態(tài)具有最小的糾纏，那么很可能典型量子熱態(tài)也具有比較小的糾纏.METTS可以通過(guò)Markov鏈Monte-Carlo采樣的方式進(jìn)行采樣，無(wú)需計(jì)算轉(zhuǎn)移到的概率是設(shè)在一步Markov采樣中，所有的態(tài)都已經(jīng)具有正確的概率那么依據(jù)上述轉(zhuǎn)移概率，在下一步迭代中，轉(zhuǎn)移到某個(gè)態(tài)的概率為

而這正是所需要的正確概率.不動(dòng)點(diǎn)的存在意味著在足夠多步之后Markov采樣將會(huì)收斂到正確的概率.另一個(gè)針對(duì)MPS的高效采樣技巧是通過(guò)逐個(gè)格點(diǎn)的方式從塌縮至

4 TD-DMRG在化學(xué)中的應(yīng)用

下面將介紹若干類TD-DMRG已經(jīng)得到成功應(yīng)用的化學(xué)科學(xué)問(wèn)題.首先將聚焦于電聲耦合模型中的非絕熱動(dòng)力學(xué)問(wèn)題，隨后將介紹基于第一性哈密頓的電子動(dòng)力學(xué)問(wèn)題.

4.1 激子及電荷轉(zhuǎn)移動(dòng)力學(xué)

激子動(dòng)力學(xué)和電荷轉(zhuǎn)移動(dòng)力學(xué)對(duì)于生物分子和新材料中的能量和電荷轉(zhuǎn)移至關(guān)重要.這類模型一般是由Frenkel-Holstein哈密頓及其變種所描述.在光合作用過(guò)程中，具有光學(xué)活性的生物分子吸收光子，形成激子.隨后，太陽(yáng)能通過(guò)分子間的激子耦合高效地轉(zhuǎn)移至反應(yīng)中心，儲(chǔ)存為化學(xué)能.從分子層面上理解這一能量轉(zhuǎn)移機(jī)制可以為優(yōu)化和設(shè)計(jì)人造光合系統(tǒng)以及光電器件提供有益的信息.2013年，Plenio等［56］使用TD-DMRG研究了典型的色素-蛋白復(fù)合物Fenna-Matthews-Olson（FMO）復(fù)合物中的激子動(dòng)力學(xué)過(guò)程.他們使用一個(gè)兩色素模型，將哈密頓變換成只具有最近鄰相互作用的等效哈密頓，然后利用TEBD算法進(jìn)行時(shí)間演化［54］.后來(lái)，Gelin和Borrelli利用TDVP-PS算法，將這一研究拓展到具有7個(gè)色素的模型以及反應(yīng)中心的電子轉(zhuǎn)移模型［50，71］.Gelin等的研究通過(guò)純化結(jié)合Bogoliubov變換的方式考慮了有限溫度效應(yīng).他們還在TD-DMRG的輔助下證明了具有靜態(tài)無(wú)序的等價(jià)色素體系的長(zhǎng)時(shí)間動(dòng)力學(xué)由振動(dòng)動(dòng)力學(xué)所控制［72］.

利用TD-DMRG也可以對(duì)有機(jī)太陽(yáng)能電池中給體/受體界面的超快激子分離問(wèn)題做類似的理論模擬［21，22］.激子分離的模型中除了激子態(tài)和電聲耦合之外一般還包括電荷轉(zhuǎn)移態(tài).Ma等［21］在不同電聲耦合強(qiáng)度α下，用TD-DMRG計(jì)算了空穴（給體部分）和電子（受體部分）密度隨時(shí)間的演化（圖4），其中電聲耦合由聲子譜密度是截止頻率）所刻畫.可見，在適中的電聲耦合強(qiáng)度下，體系中存在超快的長(zhǎng)程電荷轉(zhuǎn)移.他們將這一現(xiàn)象歸因于局域Frenkel激子態(tài)和長(zhǎng)程電荷轉(zhuǎn)移態(tài)之間的量子共振.

Fig.4 Charge density evolution of hole in the donor part and electron in the acceptor part with differentαvalues[21]

除了上述報(bào)道以外，還有針對(duì)π共軛高分子如聚對(duì)亞苯（PPP）和聚對(duì)苯撐乙烯（PPV）體系中超快激子弛豫、退相干和局域化的研究［57，73］.如果將三線態(tài)對(duì)納入到電子基組中，基于三態(tài)或更復(fù)雜的模型，TD-DMRG也可以用于研究單態(tài)裂分過(guò)程［51，74］.

這類針對(duì)激子動(dòng)力學(xué)（或電荷轉(zhuǎn)移動(dòng)力學(xué)）的理論模擬一般分為3步.首先將系統(tǒng)中沒(méi)有激子（或電荷）的零溫度純態(tài)或熱平衡態(tài)表示成一個(gè)MPS，此時(shí)系統(tǒng)中各個(gè)自由度往往沒(méi)有關(guān)聯(lián).然后，通過(guò)將產(chǎn)生算符作用在系統(tǒng)上的方式將一個(gè)激子（或電荷）注入到合適的分子上.最后就是通過(guò)合適的時(shí)間演化算法進(jìn)行系統(tǒng)的時(shí)間演化，并在每一時(shí)間步計(jì)算關(guān)鍵的物理量如電子的約化密度矩陣元.需要注意的是，根據(jù)一項(xiàng)最近使用TD-DMRG的研究，上述步驟中描述的垂直Franck-Condon激發(fā)動(dòng)力學(xué)與從振動(dòng)弛豫的激發(fā)態(tài)出發(fā)的動(dòng)力學(xué)可能完全不同［52］.此外，只要振動(dòng)自由度是由簡(jiǎn)諧振子所描述的，采用Bogoliubov變換是考慮有限溫度效應(yīng)的有效方式.

4.2 激發(fā)態(tài)動(dòng)力學(xué)及光譜

激發(fā)態(tài)動(dòng)力學(xué)與激子動(dòng)力學(xué)實(shí)際上密不可分，本文描述的激發(fā)態(tài)動(dòng)力學(xué)主要關(guān)注基于復(fù)雜勢(shì)能面的動(dòng)力學(xué)或與光譜直接相關(guān)的動(dòng)力學(xué).復(fù)雜勢(shì)能面動(dòng)力學(xué)的典型代表是紫外光激發(fā)吡嗪分子至S2態(tài)后S1/S2的內(nèi)轉(zhuǎn)換動(dòng)力學(xué).這一問(wèn)題已被包括MCTDH在內(nèi)的許多非絕熱動(dòng)力學(xué)方法所研究［75，76］.吡嗪分子共有24個(gè)正則振動(dòng)模式，電聲耦合較強(qiáng)，在S1態(tài)S2態(tài)之間存在一個(gè)錐形交叉點(diǎn).因此該體系的動(dòng)力學(xué)是量子動(dòng)力學(xué)方法的一個(gè)嚴(yán)格的基準(zhǔn)測(cè)試.基于線性響應(yīng)理論，分子的吸收光譜I(ω)可以由偶極-偶極關(guān)聯(lián)函數(shù)C(t)的傅里葉變換得到：

式中：μ為偶極算符；下角標(biāo)g代表基態(tài).在零溫度下，計(jì)算C(t)可以簡(jiǎn)化成計(jì)算時(shí)間演化的波函數(shù)Ψ(t)與初始S2態(tài)波函數(shù)Ψ(0)的重疊：

此處，Eg是基態(tài)能量，通過(guò)恰當(dāng)?shù)囟x哈密頓，可以使其為能量零點(diǎn).P&C和基于TDVP的時(shí)間演化算法都可以用于處理這一時(shí)間演化問(wèn)題［24，47，51］.這些時(shí)間演化算法與高精度的電子結(jié)構(gòu)方法導(dǎo)出的第一性參數(shù)相結(jié)合，可以模擬出與實(shí)驗(yàn)結(jié)果相吻合的光譜.

同樣的方法也可推廣到分子聚集體，如苝酰亞胺染料（PBI）的光譜［23，24］.對(duì)算法稍作調(diào)整后，也可以計(jì)算有限溫度的光譜［23］.這一調(diào)整的關(guān)鍵是在計(jì)算C(t)時(shí)，對(duì)純化后的熱平衡態(tài)同時(shí)開展兩組時(shí)間演化：

Fig.5 Calculated absorption and fluorescence spectra of PBI dimer at 298 K from finite temperature TD-DMRG[23]

4.3 電荷輸運(yùn)

共軛高分子中的電荷輸運(yùn)性質(zhì)是TD-DMRG在化學(xué)中的最早應(yīng)用之一［55，78，79］.這些工作在利用TD-DMRG描述電子運(yùn)動(dòng)的同時(shí)，利用經(jīng)典力學(xué)描述核的運(yùn)動(dòng)，形成了一種混合TD-DMRG/Ehrenfest時(shí)間演化方法.近期，我們［25］利用另一種思路即Kubo方程計(jì)算了有限溫度下有機(jī)半導(dǎo)體中的載流子電荷遷移率：

其中j是電流算符，憑借我們［80］最近開發(fā)的自動(dòng)化構(gòu)建MPO的方法，可以精確地表示成MPO［25］.圖6展示了由TD-DMRG計(jì)算得到的紅熒烯晶體遷移率和由費(fèi)米黃金規(guī)則（Fermi’s golden rule，F(xiàn)GR）和玻爾茲曼輸運(yùn)理論（Boltzmann transport theory）得到的遷移率的對(duì)比.可以看出，在跳躍極限下即轉(zhuǎn)移積分V比較小時(shí)，TD-DMRG可以回歸到FGR的解析結(jié)果.而在能帶極限下即轉(zhuǎn)移積分V比較大時(shí)，TD-DMRG的漸進(jìn)行為與Boltzmann輸運(yùn)理論相同.除此之外，TD-DMRG計(jì)算的遷移率還符合實(shí)驗(yàn)觀察到的“類能帶”傳輸行為以及負(fù)的同位素效應(yīng)［25］.

Fig.6 Carrier mobility from weak to strong electronic coupling(V)calculated by TD-DMRG,FGR(hopping limit),and Boltzmann transport theory(band limit)(A),mean free path(l mfp)of the charge carrier from weak to strong electronic coupling calculated by TD-DMRG(B)[25]

4.4 電子動(dòng)力學(xué)

盡管大多數(shù)化學(xué)中的TD-DMRG文獻(xiàn)以電聲耦合模型作為研究的對(duì)象，由于TD-DMRG是一種普適的方法，在最近的文獻(xiàn)中出現(xiàn)了純電子動(dòng)力學(xué)研究.如Chan等［38］利用改進(jìn)的TST算法和頻域DMRG算法研究了一維氫原子鏈的動(dòng)態(tài)性質(zhì)，基于第一性原理哈密頓計(jì)算了體系的態(tài)密度和復(fù)極化函數(shù)，以作為不同鍵長(zhǎng)時(shí)體系的金屬性和離域性的判據(jù).最近，F(xiàn)rahm等［39］利用TD-DMRG模擬了碘乙炔分子中基于第一性原理哈密頓的超快電荷轉(zhuǎn)移動(dòng)力學(xué)，其結(jié)果如圖7所示.初始時(shí)刻，空穴布居在碘上，隨后空穴轉(zhuǎn)移到碳上，然后又回到碘上，最后空穴再次布居在碳上，與實(shí)驗(yàn)測(cè)量結(jié)果直接符合.

Fig.7 Hole density(red)of the iodoacetylene molecule at four different points in time[39]

5 總結(jié)與展望

本文介紹了TD-DMRG的關(guān)鍵算法及近期TD-DMRG在化學(xué)問(wèn)題中的應(yīng)用.受益于現(xiàn)代DMRG精巧而簡(jiǎn)潔的MPS數(shù)學(xué)結(jié)構(gòu)，越來(lái)越多的TD-DMRG時(shí)間演化算法不斷涌現(xiàn)，這些算法各自具有不同的特點(diǎn)，需要使用者針對(duì)合適的問(wèn)題選擇合適的算法.兩種使TD-DMRG可以考慮溫度效應(yīng)的算法進(jìn)一步增強(qiáng)了TD-DMRG模擬真實(shí)體系的能力.豐富的應(yīng)用說(shuō)明，TD-DMRG可以精確地模擬許多動(dòng)態(tài)過(guò)程，如分子聚集體和生物體系中的激子動(dòng)力學(xué)、單態(tài)裂分、小分子和聚集體的一維和二維光譜、有機(jī)半導(dǎo)體的電荷輸運(yùn)以及基于第一性原理的電子動(dòng)力學(xué)等.由于矩陣乘積的數(shù)學(xué)形式簡(jiǎn)潔，未來(lái)，MPS可能會(huì)與機(jī)器學(xué)習(xí)的方法相結(jié)合用于擬合實(shí)際分子的勢(shì)能面，并在TD-DMRG的框架下發(fā)展多原子體系的量子動(dòng)力學(xué)和光譜的計(jì)算方法.另外，現(xiàn)有的TD-DMRG計(jì)算電荷傳輸遷移率的方法可以很方便地推廣到熱輸運(yùn)、自旋輸運(yùn)以及熱電轉(zhuǎn)換等問(wèn)題，具有非常廣闊的應(yīng)用空間.然而，超越一維模型是TD-DMRG模擬實(shí)際塊體材料所面對(duì)的主要挑戰(zhàn).高維模型的計(jì)算量一般比一維模型大很多，對(duì)高維模型的高效模擬仍然有賴于進(jìn)一步發(fā)展TD-DMRG方法.

TD-DMRG的一個(gè)主要不足是對(duì)最大演化時(shí)長(zhǎng)的限制.在常見的化學(xué)問(wèn)題中，體系的二分糾纏熵S隨演化時(shí)間線性增長(zhǎng)S∝t［81，82］.這意味著數(shù)值精確的模擬所需要的虛擬鍵維數(shù)隨時(shí)間指數(shù)增長(zhǎng)M∝et，或者說(shuō)如果虛擬鍵維數(shù)固定，時(shí)間演化的誤差隨時(shí)間指數(shù)增長(zhǎng)［16］.因此，總體上，TD-DMRG更適合于模擬超快過(guò)程的演化，而對(duì)長(zhǎng)時(shí)間的演化問(wèn)題仍需要發(fā)展新的算法.盡管如此，現(xiàn)有的TD-DMRG算法已經(jīng)足夠在相當(dāng)長(zhǎng)的演化時(shí)間內(nèi)給出含時(shí)薛定諤方程的可靠解.一個(gè)可以繞開這一限制的策略是通過(guò)關(guān)聯(lián)函數(shù)計(jì)算感興趣的物理量（如遷移率），因?yàn)閷?duì)于真實(shí)化學(xué)系統(tǒng)來(lái)說(shuō)關(guān)聯(lián)函數(shù)在長(zhǎng)時(shí)間應(yīng)該收斂到0.另一種策略是在頻域空間計(jì)算體系性質(zhì)［83］.在TD-DMRG中直接克服這一糾纏障礙的方法現(xiàn)在仍在發(fā)展中［84］.