胡史奇，廉超，關(guān)夢雪，孟勝,2*

1.中國科學院物理研究所北京凝聚態(tài)物理國家實驗室，北京 100190

2.量子物質(zhì)科學協(xié)同創(chuàng)新中心，北京 100190

發(fā)展動量分辨的含時密度泛函理論，實時模擬石墨烯納米帶的激發(fā)態(tài)

胡史奇1，廉超1，關(guān)夢雪1，孟勝1,2*

1.中國科學院物理研究所北京凝聚態(tài)物理國家實驗室，北京 100190

2.量子物質(zhì)科學協(xié)同創(chuàng)新中心，北京 100190

實時的密度泛函理論 (real-time Time Dependent Density Functional Theory，或稱 rt-TDDFT) 是基于含時 Kohn-Sham 方程，從實空間實時模擬體系激發(fā)態(tài)性質(zhì)的第一性原理計算方法。本文介紹如何利用自主開發(fā)的動量分辨的含時密度泛函理論計算方法和軟件 Time Dependent ab initio Package (TDAP)，研究具有較大尺寸的石墨烯納米體系的激發(fā)態(tài)性質(zhì)。以石墨烯扶手椅型納米帶 (Armchair graphenenanoribbons，AGNRs) 為例，通過直接模擬外加光場下多電子體系的波函數(shù)實時演化，得到了體系的光吸收譜以及激發(fā)態(tài)的動力學過程。這樣的計算方法大大提高了周期性體系的計算效率，且不再受微擾論的限制，實現(xiàn)了對大規(guī)模、真實體系的動力學性質(zhì)的實時模擬。

第一性原理計算；含時密度泛函；激發(fā)態(tài)模擬；石墨烯

引言

上世紀七十年代以來，伴隨著計算機技術(shù)的飛速發(fā)展，密度泛函理論[1](Density Functional Theory，DFT) 在固體物理學計算中的成功使得第一性原理計算在模擬原子尺度物理化學過程、理解預測材料的機理和性質(zhì)、發(fā)現(xiàn)設(shè)計新型功能材料等領(lǐng)域得到了廣泛的應(yīng)用。作為一種研究多電子體系電子結(jié)構(gòu)的量子力學方法，密度泛函理論利用電子體系基態(tài)總能量是電子密度的唯一泛函的假設(shè)，將復雜的電子相互作用融入到交換關(guān)聯(lián)能量項中，通過對其取局域密度或更高級的近似，得到了類似單粒子薛定諤方程形式的Kohn-Sham 方程。實際計算中，為方便求解，一般通過傅里葉變換的辦法將 Kohn-Sham 方程從實空間求解轉(zhuǎn)化到共軛動量所在的空間求解。結(jié)合自洽計算方法，通過變分求 Kohn-Sham 方程決定的總能量極小值點，得到體系在基態(tài)時的能量及性質(zhì)。密度泛函理論在凝聚態(tài)體系基態(tài)性質(zhì)中的計算取得了令人滿意的結(jié)果，是凝聚態(tài)物理、計算材料學和計算化學領(lǐng)域最常用的方法之一。

但是通常的密度泛函理論僅對于多電子體系的基態(tài)成立，對于電子激發(fā)態(tài)，原則上不存在電子密度和體系真正的多體波函數(shù)的一一對應(yīng)關(guān)系，因此無法用普通的 Kohn-Sham 方程描述。然而，E.Runge 和E.K.U.Gross 于 1984年發(fā)現(xiàn)[2]，如果初態(tài)多體波函數(shù)給定，那么隨時間變化的電子密度和隨時間變化的體系多體波函數(shù)也具有嚴格的一一對應(yīng)關(guān)系，此即含時密度泛函理論。依據(jù)含時密度泛函理論，體系的激發(fā)態(tài)性質(zhì)原則上可以嚴格求解，在實際應(yīng)用中則需要對時間相關(guān)的交換關(guān)聯(lián)泛函取近似，并發(fā)展實際可行的算法。

一種實際辦法就是實時求解含時密度泛函方程(real-time Time Dependent Density Functional Theory，或稱 rt-TDDFT)，直接模擬體系的激發(fā)態(tài)性質(zhì)。相比于傳統(tǒng)的基于準粒子修正描述激發(fā)態(tài)行為的格林函數(shù)方法，rt-TDDFT 包含了所有的非線性效應(yīng)，擺脫了微擾論要求的必須在平衡態(tài)附近應(yīng)用的限制。同時，rt-TDDFT 方法提供了電子波函數(shù)實時的演化過程，結(jié)合離子動力學，能夠直觀地給出離子-電子多體體系隨時間的演化路徑。因此，rt-TDDFT 成為計算模擬強場物理、超快物理過程的有效工具。

然而，受限于計算效率，傳統(tǒng)的 rt-TDDFT 大多只能在動量空間的單點上進行計算，無法實現(xiàn)體系的動量分辨。這一限制意味著傳統(tǒng)的 rt-TDDFT 無法對周期性的體系進行有效大規(guī)模模擬，只能限于原子團簇等零維非周期體系。同時，傳統(tǒng)的 rt-TDDFT 大都只考慮了電子的演化進程，而不涉及離子動力學對電子狀態(tài)的影響。這些都大大限制了 rt-TDDFT 理論的應(yīng)用范圍，從而阻礙了 rt-TDDFT 在物理、材料計算領(lǐng)域的應(yīng)用和發(fā)展?；谏鲜霈F(xiàn)狀，我們發(fā)展了一套具有動量空間分辨能力，并且考慮離子動力學過程的 rt-TDDFT 計算方法和軟件 TDAP (Time Dependent ab initio Package)[3-4]。基于 TDAP，我們不僅實現(xiàn)了周期性體系的激發(fā)態(tài)性質(zhì)的計算模擬，同時計算效率也大大提高。我們選取了近年來的熱點體系—扶手椅型石墨烯納米帶 (Armchair graphenenanoribbons 或稱AGNRs) 作為研究對象，通過外加光場的方式，模擬了不同尺寸 AGNRs 的光吸收行為以及不同吸收模式下的激發(fā)態(tài)動力學過程。我們發(fā)展的模擬方法實現(xiàn)了電子激發(fā)模式在動量空間的分辨，擺脫了傳統(tǒng)微擾論的限制，同時也提高了計算效率。

1 原理方法

1.1 含時 Kohn-Sham 方程

多電子體系的哈密頓量為:

其中第一、二項代表了電子與離子的運動動能，第三項、第四項分別代表了電子間、離子間的庫倫相互作用，第五項代表了電子-離子之間的相互作用，最后一項則描述了外勢場。此時電子體系的量子行為一般采用薛定諤方程進行描述，在 rt-TDDFT 理論中，多體電子體系的薛定諤方程可以在一定程度上簡化為單粒子形式的含時 Kohn-Sham (Time-dependent Kohn-Sham，TDKS) 方程[3]：

其中，φj(r,t) 為第j個電子的單粒子波函數(shù)，第二項描述了外場的貢獻，第三項表示了電子的庫倫相互勢能項，第四項描述了離子提供的勢場，而電子-電子之間的交換關(guān)聯(lián)作用，被寫入最后一項交換關(guān)聯(lián)能中。

1.2 電子波函數(shù)演化求解

對于體系激發(fā)態(tài)性質(zhì)的描述可以通過求解 TDKS方程得到，各時刻的電子密度矩陣由單粒子波函數(shù)表示。同時，t+dt時刻單電子波函數(shù)可以通過t時刻的波函數(shù)演化得到[4]：

其中Sk(t') 為動量k標記的波函數(shù)的交疊矩陣，Hk(t')為動量k標記的t' 時刻的哈密頓量，二者又依賴于此時的電子密度。因此，可以用自洽計算描述電子體系的演化。演化過程中，電子各時刻的不同k點能級變化可以通過對該時刻哈密頓量進行對角化求解本征值得到[4]：

同時，躍遷占據(jù)數(shù)也由演化的 TDKS 波函數(shù)和各時刻哈密度量對角化的本征態(tài)基矢的內(nèi)積表示：

2 離子動力學過程

根據(jù) Ehrenfest 理論，在忽略離子間的交換關(guān)聯(lián)相互作用后，離子的運動狀態(tài)可以利用經(jīng)典的牛頓方程來描述[3]：

3 原子軌道基矢

圖1 TDAP 計算框架Fig.1 Flowchart of TDAP algorithm

TDAP 是基于 Siesta 計算軟件。在計算中，我們采用數(shù)值原子軌道基矢 (Numerical Atomic orbital 或稱 NAOs) 作為展開。相比于平面波基矢，原子軌道基矢由于其局域特性大大提高了收斂速度。因此，原子軌道基矢的使用在保證精度的同時提高了計算的速度。對于有Nn個 TDKS 軌道，No個原子軌道基矢的體系，計算速率呈現(xiàn)NnO(No2) 的依賴關(guān)系，對激發(fā)態(tài)的模擬時長也得到提高。

4 激發(fā)態(tài)模擬

計算中，對體系的激勵來自于模擬的外界光場，這樣的外加光場在程序上的實現(xiàn)是通過對哈密頓量中的外勢進行改造，添加外界光場產(chǎn)生的勢。這里，主要考慮光的電場分量[4]：

對于周期性體系而言，周期性電荷密度分布的限制導致電場的施加必須以存在真空層為前提，因此計算中，電場方向應(yīng)該沿體系的非周期性方向。實際模擬中，主要采用兩種外加電場形式，對材料吸收特性的模擬，采用包含各個頻率分量的階躍場E(t)=E(1-Θ(t))；對激發(fā)態(tài)電子動力學過程的考察采用特定頻率的正弦場：E(t)=Esin (2πft+φ)。根據(jù)電偶極矩與外加電場的關(guān)系：P(ω)=ε(ω)E(ω)，施加電場激勵后，通過探測體系不同時刻電偶極矩的響應(yīng)可以得到體系的介電性質(zhì)：

進而可以得到體系吸收譜性質(zhì)：Abs(ω)∝ Im(ε(ω))，以及能量損失譜 (EELS) 性質(zhì)：EELS(ω)∝-Im(1/ε(ω))。

5 結(jié)果與討論

圖2 石墨烯納米帶光場激發(fā)示意圖Fig.2 The schematic showing the GNRs under an electric field

首先計算石墨烯納米帶的基態(tài)電子結(jié)構(gòu)。我們采用 Siesta 軟件包，采用模守恒贗勢和局域密度近似，在數(shù)值原子軌道基的基礎(chǔ)上對石墨烯納米帶的結(jié)構(gòu)進行了弛豫。其中k點的選取為 1×1×24，能量收斂精度為 10-4eV。在此基礎(chǔ)上，我們利用 TDAP，通過模擬沿寬度方向的外界光場，計算了石墨烯納米帶的吸收譜、能量損失譜等性質(zhì)。在對吸收譜與能量損失譜的模擬中，我們采用階躍電場，場強大小為0.25V/?，模擬時間為 24fs，同時選取 1×1×32 的動量空間分布。隨后，我們將 TDAP 的模擬結(jié)果與 DFT、TDDFT 線性響應(yīng)的結(jié)果進行了對比，測試了程序的可靠性。在得到吸收性質(zhì)的基礎(chǔ)上，對激發(fā)態(tài)電子的動力學過程進行了探討。得到了不同吸收模式下的激發(fā)動力學。對于激發(fā)動力學的模擬，我們采用正弦電場，場強大小為 0.125V/?，模擬時間為 50fs，選取1×1×32 的動量空間分布。

5.1 石墨烯納米帶的光吸收譜與能量損失譜

圖3給出了不同寬度不同能量區(qū)域石墨烯納米帶的吸收譜特性，根據(jù)不同能量區(qū)域，光吸收存在兩個模式：可見光區(qū)域的模式和高能量區(qū)域的模式。兩類模式的共振吸收能量均隨著寬度的增加而紅移。其中，高能量區(qū)域的吸收模式為主導，其共振吸收能量大約在 5-7eV。對比石墨烯納米帶的計算結(jié)果與DFT、DFT+GW 計算得到的石墨烯吸收譜[6,7]，可以發(fā)現(xiàn)，隨著寬度增加，兩類吸收模式以不同的方式趨向于無限大尺寸的二維石墨烯材料的行為。其中，對可見光區(qū)域，隨寬度增加，各吸收峰逐漸展寬最終形成趨向于石墨烯 2.3% 的吸收平臺區(qū)。對于高能量區(qū)域模式，隨著寬度的增加，其吸收峰能量逐漸趨向于石墨烯的π能帶躍遷。相較于基態(tài)的 DFT 計算，TDAP 的計算更趨向于包含準粒子修正的 DFT+GW的計算結(jié)果，從側(cè)面也印證了 TDAP 對激發(fā)態(tài)特性描述的合理性。

圖4中對比了對應(yīng)不同波矢的各寬度石墨烯納米帶電子能量損失譜的 TDAP 計算結(jié)果與利用頻率空間 TDDFT 線性響應(yīng)理論得到的石墨烯電子能量損失譜[5]。可以看出，不同寬度的石墨烯納米帶中的激發(fā)態(tài)，在一定意義上可以看做是石墨烯對應(yīng)不同波失的激發(fā)模式。二者在不同能量區(qū)域的分布均呈現(xiàn)出一致性。其中 5-7eV 的區(qū)域?qū)?yīng)于石墨烯的π能帶帶間等離激元。二者分布趨勢的一致性不僅是石墨烯納米帶與無限大石墨烯片的電子結(jié)構(gòu)密切關(guān)聯(lián)的體現(xiàn)，也說明了 TDAP 程序在激發(fā)態(tài)計算中的準確性。同時，體系的區(qū)別導致二維石墨烯內(nèi)部電子與石墨烯納米帶帶間庫倫相互作用存在差異，也使得二者在能量損失譜中細節(jié)上的不同。比如可見光區(qū)域，石墨烯納米帶的能量損失譜呈現(xiàn)一個峰值分布的振蕩，而石墨烯的能量損失譜呈現(xiàn)均一的的平臺。

圖3 不同寬度石墨烯納米帶與石墨烯的吸收譜，其中 N 為整個石墨烯納米帶沿寬度方向碳原子二聚體的數(shù)目Fig.3 Absorption spectrum of grapheme and AGNRs with different widths,N is the number of dimer lines along the ribbon width.

圖4 (a) 石墨烯與 (b) 石墨烯納米帶的電子能量損失譜Fig.4 Electron Energy Loss Spectroscopy (EELS) of (a)grapheme and (b) AGNRs.

5.2 石墨烯納米帶激發(fā)態(tài)電子動力學

吸收譜和電子能量損失譜對TDAP與DFT+GW以及TDDFT線性響應(yīng)理論計算的比較驗證了自主發(fā)展的程序在計算激發(fā)態(tài)時結(jié)果的合理性。在此基礎(chǔ)上，有別于基態(tài)的計算，借助TDAP計算效率的優(yōu)勢，我們進一步考察了石墨烯納米帶體系不同吸收模式下的激發(fā)態(tài)電子動力學過程。動力學過程能夠更為直觀的給出激發(fā)態(tài)電子波函數(shù)隨時間變化的行為。

圖5 碳原子二聚體數(shù)目為 6 的石墨烯納米帶 (a) 可見光吸收模式與 (b) 高能吸收模式激發(fā)電子總量 (藍點)，體系總能量 (紅點) 與電偶極矩 (綠線) 隨時間變化，其中黑線為外加電場隨時間的變化。Fig.5 The total excited electrons (blue),total energy (red) and dipole (green) for (a) visible light mode and (b) high energy mode on6-AGNRs,black line shows the electric field.

圖5 給出碳原子二聚體數(shù)目N為6的石墨烯納米帶(6-AGNR) 可見光區(qū)域與高能量區(qū)域不同吸收模式對應(yīng)的激發(fā)量，總能量以及電偶極矩隨時間的變化。其中黑線描述了外界光場的振蕩。可以看出，對于不同的吸收模式，電子激發(fā)量、總能量以及電偶極矩隨時間變化有所差別，意味著不同的激發(fā)態(tài)電子行為。

圖6 不同時刻下，碳原子二聚體數(shù)目為6的石墨烯納米帶的(a)可見光吸收模式與 (b) 高能吸收模式的電子激發(fā)量在不同k點上的分布。其中紅色表示占據(jù)數(shù)增加，藍色表示占據(jù)數(shù)減少，圓圈大小代表了變化數(shù)的多少。Fig.6 Theexcited electrons with k resolution at different times for (a) visible light mode and (b) high energy mode on 6-AGNRs.Red (blue) represents the increase (decrease) of electrons andcircle size is proportional to the population.

借助 TDAP 對不同動量的分辨能力，圖6展示了碳原子二聚體數(shù)目為6的石墨烯納米帶不同吸收模式激發(fā)過程中，能帶上不同k點電子占據(jù)數(shù)隨時間的變化，對比兩類的吸收模式，電子的激發(fā)躍遷在k空間分布存在差別。對于可見光模式，電子的激發(fā)趨向于個體k點，相比之下，高能吸收模式中，電子的激發(fā)分布于多個k點，趨向于集體性質(zhì)的k點分布。因此，可見光吸收模式有著個體激發(fā)的特征，而高能模式包含更多集體激發(fā)的特點。

6 結(jié)論

本文從計算原理與算法的角度，首先介紹了含時密度泛函理論在計算體系激發(fā)態(tài)性質(zhì)方面的應(yīng)用。作為基于 TDKS 演化方程的計算理論，rt-TDDFT 包含了所有的非線性效應(yīng)，提供了電子波函數(shù)隨時間變化的行為，能夠直觀地給出體系隨時間的演化路徑。同時，受限于計算效率，傳統(tǒng)的 rt-TDDFT 大都無法實現(xiàn)體系的動量分辨，并且無法考慮離子運動對電子態(tài)變化的影響。這樣的限制意味著其只能限于原子團簇等零維非周期體系的電子演化計算?；诖耍覀儼l(fā)展了包含動量分辨，同時考慮離子動力學過程的 rt-TDDFT 計算程序 TDAP。我們利用 TDAP 計算了石墨烯納米帶的吸收特性，以及不同吸收模式下的電子動力學過程。通過與已有體的計算結(jié)果進行對比，計算結(jié)果不僅驗證了程序的可靠性與高效性，同時也給出了石墨烯納米帶體系激發(fā)態(tài)的動力學圖像。本文能夠?qū)ぐl(fā)態(tài)的計算模擬、以及材料激發(fā)特性的理解與應(yīng)用提供參考價值。

[1]P.Hohenberg and W.Kohn,Phys.Rev.136,B864 (1964).

[2]E.Runge and E.K.U.Gross,Phys.Rev.Lett.52,997(1984).

[3]S.Meng and E.Kaxiras,J.Chem.Phys.129,054110(2008).

[4]C.Lian,S.-Q.Hu,M.-X.Guan,and S.Meng,arXiv:1702.08163v1(2017).

[5]C.Vacacela Gomez,M.Pisarra,M.Gravina,J.M.Pitarke,and A.Sindona,Phys.Rev.Lett.117,116801(2016).

[6]O.V Sedelnikova,L.G.Bulusheva,and a V Okotrub,J.Chem.Phys.134,244707 (2011).

[7]P.E.Trevisanutto,M.Holzmann,M.C?té,and V.Olevano,Phys.Rev.B.81,121405(2010).

Momentum-Resolved rt-TDDFT Algorithm for Modeling Electronic Excitation in Graphene Nanoribbons

Hu Shiqi1,Lian Chao1,Guan Mengxue1,Meng Sheng1,2*

1.Beijing National Laboratory for Condensed Matter Physics,Insitute of Physics,Chinese Academy of Sciences,Beijing 10018,China
2.Collaborative Innovation Center of Quantum Material,Beijing 100871,China

Real-time (rt) time dependent density functional theory (TDDFT) is an efficient ab initio method to study electronic dynamics in real time and real electron-nuclear space.The method is based on time dependent evolution of Kohn-Sham (TDKS) orbitals.Here we report first-principles simulations on electronic excitation dynamics in grapheme nanoribbons (AGNRs) using our home-built Time Dependent ab initio Package (TDAP).We obtain optical absorption spectrum and electron energy loss spectrum (EELS) via applying an external field to AGNRs and then analyze the dynamic evolution of electronic excited states.This algorithm represents an effective way to simulate electronic excitations in periodic systems beyond the limit of perturbation theory.

ab initio method; TDDFT; electronic excitation; graphene

10.11871/j.issn.1674-9480.2017.02.007

國家重點研發(fā)計劃 (2016YFA0300902、2015CB921001)；國家自然科學基金 (11774396)

*通訊作者：孟勝(smeng@iphy.ac.cn)

2017年1月13日

胡史奇：中國科學院物理研究所，博士，主要研究方向為低維材料的激發(fā)態(tài)電子性質(zhì)模擬。

E-mail:isqhu@iphy.ac.cn

廉 超：中國科學院物理研究所，博士，主要研究方向為激發(fā)態(tài)的電子動力學行為模擬。

E-mail:chaolian@iphy.ac.cn

關(guān)夢雪：中國科學院物理研究所，博士，主要研究方向為低維材料的激發(fā)態(tài)電子性質(zhì)模擬。

E-mail:mxguan@iphy.ac.cn

孟 勝：中國科學院物理研究所，量子物質(zhì)科學協(xié)同創(chuàng)新中心，研究員，博士生導師，主要研究方向為激發(fā)態(tài)的電子動力學行為，新能源器件物理以及表面水結(jié)構(gòu)研究。

E-mail:smeng@iphy.ac.cn