呂程燁 陳英煒2) 謝牧廷 李雪陽2) 于宏宇2) 鐘陽2) 向紅軍2)3)?

1)(復(fù)旦大學(xué)物理學(xué)系與計算物質(zhì)科學(xué)研究所,計算物質(zhì)科學(xué)教育部重點實驗室,上海 200433)

2)(上海期智研究院,上海 200030)

3)(人工微結(jié)構(gòu)科學(xué)與技術(shù)協(xié)同創(chuàng)新中心,南京 210093)

1 引言

“第一性原理計算”是一種基于量子力學(xué)的計算方法,通常以密度泛函理論(DFT)為基礎(chǔ),用于對材料和物質(zhì)的性質(zhì)進(jìn)行數(shù)值計算.該方法的優(yōu)點在于: 原則上不需依賴任何經(jīng)驗參數(shù)或特定的模型假設(shè),只需要知道構(gòu)成物質(zhì)的原子的種類和排布就能從“第一性原理”出發(fā)得到計算結(jié)果.因此,第一性原理計算方法在物理學(xué)、化學(xué)、材料科學(xué)、地質(zhì)科學(xué)等領(lǐng)域中有著廣泛的應(yīng)用.

DFT 最早由Hohenberg 和Kohn[1,2]在1964年提出,他們證明了如下兩條定理(即HK 定理): 基態(tài)非簡并的多體電子系統(tǒng)的性質(zhì)由基態(tài)電子密度唯一決定;且系統(tǒng)能量可視為密度的泛函,當(dāng)密度取基態(tài)密度時,能量取最小值即基態(tài)能.然而,這兩條定理并沒有給出能量泛函的具體形式,所以DFT 在此時仍然只是一個形式理論.為使DFT 成為一個可實際使用的理論,次年即1965年,Kohn和Sham[3]提出了著名的Kohn-Sham (KS)方法:他們引入了一個虛擬的無相互作用電子系統(tǒng)(即KS 系統(tǒng)),該系統(tǒng)的基態(tài)電子密度與真實系統(tǒng)相同,因此可以通過求解KS 系統(tǒng)來得到真實系統(tǒng)的基態(tài)電子密度和基態(tài)能量,從而將復(fù)雜的多電子問題轉(zhuǎn)化為可解的單電子問題.在這個方法中,每個電子都遵循相同形式的哈密頓量VKS(r),其中勢能分為三項,VKS(r)=Vext(r)+VH(r)+Vxc(r),分別為外勢(包括原子核對電子的庫侖吸引勢Vn)、電子哈特里(Hartree)勢和交換關(guān)聯(lián)勢;交換關(guān)聯(lián)勢的具體形式未知,計算中需使用合理的交換關(guān)聯(lián)泛函近似;其中Hartree 勢和交換關(guān)聯(lián)勢都依賴于密度,因此對KS 體系的求解需使用自洽場方法[3].由于密度泛函理論在后世的應(yīng)用幾乎都基于KS 方法,所以二者又常被合稱為KS-DFT.

通過施加電場或磁場來調(diào)控物質(zhì)性質(zhì)一向是科學(xué)研究的焦點,如在鐵電學(xué)[4]、鐵磁學(xué)[5]、磁光學(xué)[6]、自旋電子學(xué)[7]、谷電子學(xué)[8]、超導(dǎo)材料[9]、多鐵材料[10]、低維材料[11]、巨磁電阻效應(yīng)[12]、量子霍爾效應(yīng)[13]、聲子霍爾效應(yīng)[14]等領(lǐng)域中都需要外場的作用.為了理解和預(yù)言外場下的物性,亟需發(fā)展外場下周期性體系的第一性原理計算方法.然而,HK 定理的證明嚴(yán)格依賴于多體哈密頓量的形式和基態(tài)的存在性,因此DFT 并不能直接應(yīng)用于外加勻強(qiáng)電磁場的情形,這無疑是一個重大的挑戰(zhàn).傳統(tǒng)的解決辦法是將外場視為微擾,利用基于線性響應(yīng)的格林函數(shù)方法[15]或密度泛函微擾理論(DFPT)[16,17]進(jìn)行計算,但這些方法不能適用于有限場(包括強(qiáng)場)的情形.

本文將著重回顧近幾十年來關(guān)于有限外場下第一性原理計算方法的進(jìn)展和應(yīng)用.若無特殊說明,本文公式均選用原子單位制,即?=e=me=1 (其中 ? 為約化普朗克常數(shù),e為電子電荷量,me為電子質(zhì)量),且本文默認(rèn)采用玻恩-奧本海默近似.

2 有限電場下的DFT 計算

在經(jīng)典電磁理論中,一般用極化強(qiáng)度的變化來描述絕緣體(或稱電介質(zhì))對外電場的響應(yīng).同樣地,在后續(xù)部分,我們將看到有限電場下的第一性原理計算也直接依賴于對極化強(qiáng)度的計算.盡管外電場產(chǎn)生的外電勢不符合無限大晶體的周期性邊界條件,但極化強(qiáng)度卻可以在原胞中明確定義.因此,在討論有限電場下的第一性原理計算時,必須首先討論極化強(qiáng)度的計算.值得注意的是,本文討論的極化強(qiáng)度以及對外電場響應(yīng)的計算只適用于絕緣體.對于周期性的金屬體系,目前尚無可靠的第一性原理方法.

2.1 現(xiàn)代電極化理論

根據(jù)經(jīng)典理論,當(dāng)晶體極化時,晶體分子的正負(fù)電荷中心會發(fā)生偏移,從而形成電偶極,極化強(qiáng)度矢量就是單位體積內(nèi)的電偶極矩.但以上定義只是宏觀層面的,而微觀層面的電極化定義和計算直到20 世紀(jì)90年代才由Resta 等[18–22]發(fā)展成熟,這一理論后來被稱為“現(xiàn)代電極化理論(MTP)”.對于我們關(guān)心的無限大晶體,需要尋找的是極化強(qiáng)度的“體”定義,即不依賴于晶體表面的極化性質(zhì).一種直觀的定義是把極化強(qiáng)度定義為原胞內(nèi)的平均電偶極矩(其中?是原胞體積,ρ(r)為電荷密度),但這勢必會依賴于原胞位置和形狀的選取,從而導(dǎo)致結(jié)果不唯一,如圖1所示.況且即使從麥克斯韋方程組出發(fā),?·P(r)=-ρ(r),電荷密度ρ(r)也不能唯一決定極化強(qiáng)度分布P(r),因為P(r)加上任意一個無散場都 滿足此方程.以上種種說明周期性體系的極化強(qiáng)度不能由電荷密度唯一決定.

圖1 兩種原子構(gòu)成的系統(tǒng)的不同原胞選擇[21],兩原子電荷分別為 Z1=+e (空心圓)和 Z2=+3e (陰影圓)(a),(b)包含了完整的原子,但相對位置不同;(c)原胞由一個完整的+e 電荷和4 個+3e/4 電 荷組成Fig.1.Possible choices of unit cell for a system composed of two types of atoms having ionic charges Z1=+e (open circles)and Z2=+3e (shaded circles)[21]: (a),(b)Unit cell is specified by two complete basis ions,but in different relative orientations;(c)unit cell is specified by“split basis”consisting of one complete+e charge and four charges+3e/4 in a symmetric arrangement.

MTP 的高明之處在于: 它并不直接定義電極化強(qiáng)度本身,而是定義極化強(qiáng)度在晶體絕熱演化下初末態(tài)之間的變化量[18].事實上實驗學(xué)家測量的也正是極化強(qiáng)度的變化量而非極化強(qiáng)度本身,如鐵電材料的自發(fā)極化是通過測量電滯回線得到的.極化強(qiáng)度的變化量可以被定義為

其中 i和f 分別指代初、末態(tài),總極化強(qiáng)度P可認(rèn)為是離子極化強(qiáng)度Pion和電子極化強(qiáng)度Pel的和,離子極化強(qiáng)度定義為

其中Zμ和uμ分別為離子μ的電荷數(shù)和位矢.

在絕熱近似下,對應(yīng)本征能Enk的電子本征態(tài)仍滿足布洛赫函數(shù)形式,即ψnk(r)=eik·runk(r),其中unk(r)=unk(r+R)是布洛赫函數(shù)的周期部分,n為能帶指標(biāo),k取第一布里淵區(qū)中的點,R為晶格的正格矢.由一階絕熱微擾論給出

其中對n求和即對所有 共M條占據(jù)帶求和,f是自旋簡并度,負(fù)號來自于電子帶負(fù)電.由此可以給出與演化路徑無關(guān)的極化強(qiáng)度:

其中An(k)=i〈unk|?kunk〉 就是k空間上的Berry聯(lián)絡(luò).也常將(3)式寫成Berry 相位的形式:

其中ai是原胞基矢,

是沿倒格矢基矢bi方向的Berry 相位,是第一布里淵區(qū)的體積.極化強(qiáng)度正比于Berry 相位的事實也暗示了其應(yīng)當(dāng)是多值的,因為Berry 相位在相差 2π 的意義上是規(guī)范不變的,即2πm,m為任意整數(shù).所以并不能直接比較(2)式和(3)式給出的極化強(qiáng)度的大小,但是兩個狀態(tài)的極化差可以通過構(gòu)造一個絕熱演化路徑唯一確定.此外(3)式還有一個更具物理直觀的詮釋: 定義Wannier 中心

其中Wannier 函數(shù)

需要指出,前面討論的MTP 是一個單電子理論,真實體系的極化強(qiáng)度應(yīng)由多體理論給出[25,26].Resta 等[26,27]通過巧妙構(gòu)造與周期性邊界條件適配的多體位置算符給出了多體極化強(qiáng)度的正確定義,其在多體波函數(shù)取Slater 行列式時退化到單體的MTP.

2.2 固定電場方法和電熱焓

現(xiàn)在需要研究電子在有限大外電場E中的性質(zhì).考慮KS-DFT 框架,晶體中電子的哈密頓量寫作[28]:

當(dāng)存在外電場時,我們認(rèn)為體系“基態(tài)”仍處于極化的布洛赫態(tài)ψnk(r)=eik·runk(r),此時體系的總能量可以用電熱焓F來描述,體系“基態(tài)”即對應(yīng)于取電熱焓為極小值時的波函數(shù).考慮體系對外電場的響應(yīng),電熱焓F可以定義為[36]

其中P=Pion+Pel是總極化強(qiáng)度,離子極化強(qiáng)度Pion和電子極化強(qiáng)度Pel分別由(2)式和(3)式給出.需要指出的是,ψnk并不是哈密頓量的本征態(tài),但可以被視為對單粒子密度矩陣的表示,且該密度矩陣仍保留有周期性[32].在對(6)式進(jìn)行最優(yōu)化時,可認(rèn)為外電場E是一個常量,因而該算法又被稱作“固定電場方法”.物理上這對應(yīng)于給材料施加一個恒定的電壓(暫不考慮形變),即將材料兩端分別與電源的正負(fù)極相連形成閉路.基于電熱焓需要取到最小值,數(shù)值上可以采用預(yù)條件處理的共軛梯度法來求解該最優(yōu)化問題[35].

我們有時需要對原子的位置以及原胞的大小進(jìn)行更新.這里以優(yōu)化原子位置為例: 借助于Hellmann-Feynman 定理,可以計算離子受力Fμ=,該式第1項與傳統(tǒng)DFT計算一致,第2 項則是外電場對離子的作用力,其中Zμ為μ離子的電荷數(shù).另一種常見的離子受力算法需利用玻恩有效電荷,其中μ原子的玻恩有效電荷張量定義為,一般通過極化強(qiáng)度關(guān)于μ原子沿β方向位移uμ,β的有限差分得到.基于受力,可以通過共軛梯度法或其他優(yōu)化方法來優(yōu)化原子位置,使總能最小.值得一提的是,我們可以通過犧牲一定精度來提高效率.Fu 和Bellaiche[36]直接利用零場下的玻恩有效電荷來計算加電場后的原子受力,以結(jié)構(gòu)優(yōu)化,而舍棄了對電子波函數(shù)的優(yōu)化,也可以計算出體系的由晶格主導(dǎo)的壓電和介電響應(yīng)系數(shù),且與實驗結(jié)果符合得較好;同時,利用這種方法還可以計算體系的電滯回線[37,38];此外,利用零場玻恩有效電荷來近似有限場下玻恩有效電荷,在后文將要介紹的有效哈密頓量方法中也十分普遍.

但該算法也存在一系列問題.如果k空間取點過密,則F將不存在極小值,從而導(dǎo)致算法失效;而在一般的DFT 計算中,k空間取點越密意味著計算越精確.這個缺陷的背后是兩種特征長度的比較: 其一是Zener隧穿距離Lt=Egap/|E|,其 中Egap為能隙大小;其二是Lp=2π/|?k|,其中?k為k點采樣間距,該特征長度可以被視為施加周期性邊界條件的長度.為能正確求解“基態(tài)”波函數(shù),我們需要求Lp

利用固定電場方法,可以對體系的介電性、壓電性、多鐵性等許多相關(guān)的物理量進(jìn)行計算模擬.比如,玻恩有效電荷也可以通過改變電場大小用對受力的有限差分來計算.固定電場方法結(jié)合有限差分還可以計算極化率,以及介電常數(shù)εαβ=δαβ+χαβ,其中ε0是真空介電常數(shù).如果在計算中使離子固定不動,將得到電子貢獻(xiàn)ε∞;如果同時允許離子和電子的弛豫,則將得到靜態(tài)介電常數(shù)εstatic.同樣可以將該方法推廣至二階極化率計算中應(yīng)固定離子位置.利用不同贗勢方法對一些實際體系的介電性的計算結(jié)果可見表1,與實驗符合得較好[42].當(dāng)然,這些零場下的介電性相關(guān)的物理量,也可以通過DFPT 進(jìn)行計算[43].

表1 用固定電場方法計算的一些III-V 半導(dǎo)體介電性質(zhì)與實驗的比較[42],其中玻恩有效電荷張量在材料對稱性下退化為標(biāo)量,且 d123 定義為 /2,LDA 和PBE 是計算時使用的交換關(guān)聯(lián)泛函近似Table 1.Computed dieletric properties of some III–V semiconductors by means of fixed-E method compared to experiment[42],Born effective charge tensor degenerates to a scalar due to the symmetry of the material and d123 is defined as /2,LDA and PBE are different exchange-correlation functional approximations used in calculation.

表1 用固定電場方法計算的一些III-V 半導(dǎo)體介電性質(zhì)與實驗的比較[42],其中玻恩有效電荷張量在材料對稱性下退化為標(biāo)量,且 d123 定義為 /2,LDA 和PBE 是計算時使用的交換關(guān)聯(lián)泛函近似Table 1.Computed dieletric properties of some III–V semiconductors by means of fixed-E method compared to experiment[42],Born effective charge tensor degenerates to a scalar due to the symmetry of the material and d123 is defined as /2,LDA and PBE are different exchange-correlation functional approximations used in calculation.

壓電效應(yīng)是指當(dāng)體系施加電場后,會產(chǎn)生原胞的形變.這一效應(yīng)具有廣泛場景,比如壓電濾波器、加速器和壓力傳感器等.可以用壓電張量來描述壓電效應(yīng),其可定義為應(yīng)力張量對場強(qiáng)的導(dǎo)數(shù)eαβγ=,其中σαβ=為應(yīng)力張量,η為應(yīng)變張量,Vi=E·ai為電壓.但當(dāng)考慮體系存在形變時,即ai→(1+η)ai時,固定電壓Vi和固定電場E不再等同,見圖2.由于實驗上一般是在材料兩側(cè)施加固定電壓,因此需要基于固定電壓法來進(jìn)行計算模擬.此時電場強(qiáng)度會與應(yīng)變張量耦合,可以定義約化電場強(qiáng)度E′=E(1+η)和約化極化強(qiáng)度P′=(1+η)-1P使之解耦[25,26],這樣定義的電場強(qiáng)度才是考慮形變時的電熱焓的自然變量:

圖2 固定電壓法和固定電場法的差異[45].當(dāng)材料發(fā)生應(yīng)變 η時,如果保持電壓?V不變,電場從E=?V/d 變 化為 E=?V/[(1+η)d] ;如果保持電場 E 不變,電壓從 ?V 變?yōu)?1+η)?VFig.2.Differences between fixed-E method and fixedvoltage method[45].When a strain η is applied to the material,electric field will change from E=?V/d to E=?V/[(1+η)d]if voltage is held fixed,or voltage will change from ?V to (1+η)?V if electric field is held fixed.

其中?(η)是形變后的原胞體積.通常把基于固定電壓得到的壓電張量稱為恰當(dāng)?shù)?proper)壓電張量,而基于固定電場得到的為不恰當(dāng)?shù)?improper)壓電張量,Vanderbilt[44]曾推導(dǎo)過兩者之間的轉(zhuǎn)換關(guān)系,得到.已經(jīng)有許多工作運用上述方法,來計算模擬不同體系的壓電系數(shù)[32,45].注意到,壓電系數(shù)還可以表述為eαβγ=?Pα/?ηβγ,故也可以通過DFPT 進(jìn)行計算[43].

固定電場方法同樣可以用來研究材料的多鐵性.所謂“多鐵性”,是指一種材料同時具備鐵電和鐵磁的性質(zhì),該類材料被視為新型多功能磁電器件和高性能信息存儲與處理器件的理想候選材料[46].在多鐵性物質(zhì)的研究中,主要關(guān)注如何通過電場調(diào)控磁性和通過磁場調(diào)控鐵電性,所以通常用磁電張量來描述體系的多鐵性.借助固定電場方法,可以利用有限差分來計算磁電張量.例如,Malashevich 等[47]成功利用了固定電場方法計算了多鐵材料Cr2O3中磁矩關(guān)于電場的線性響應(yīng),從而得到了磁電張量.

最后需要強(qiáng)調(diào)的是,盡管該方法計算得出的結(jié)果與實驗符合良好,但基于一般KS-DFT 的方法在原則上并不能給出真實的極化強(qiáng)度,這點最早由Gonze 等[48]注意到.該論斷的最直接證據(jù)是KS-DFT 本就不能給出真實的多體波函數(shù),而只能給出正確的基態(tài)密度,但我們知道光憑電荷密度是無法確定極化強(qiáng)度的[49].而根據(jù)Gonze 等[48]提出的“密度-極化泛函理論”,倘若想同時確定密度和極化強(qiáng)度,則必須引入有效電場EKS,但它并不等于真實的外電場E,二者的差值EKS-E被定 義為“交換關(guān)聯(lián)電場”.

2.3 固定電位移方法和內(nèi)能泛函

在實驗研究中,除了將材料接入閉路即對應(yīng)固定電場方法外,還可以讓其處于開路,這可以通過使測量電路的電阻遠(yuǎn)大于材料樣品的電阻近似實現(xiàn).在這種狀態(tài)下,材料兩端的自由電荷不會發(fā)生轉(zhuǎn)移,從而使得電位移矢量保持恒定.基于這種考慮,Stengel 等[50]提出了“固定電位移方法”,其核心是最小化如下泛函(取高斯單位制):該泛函事實上就是外加電場時原胞的內(nèi)能.Hong和Vanderbilt[51]以PbTiO3為例,測試了電場隨電位移矢量的變化,電位移矢量隨極化強(qiáng)度的變化以及極化強(qiáng)度隨電場的變化,如圖3 所示,可見電位移矢量與極化強(qiáng)度之間是一一映射的.

圖3 施加場的方向約束在[001],[110] 或[111] 方向時,PbTiO3 中形式為 ε (D)(ε 為外電場)(a)、D (P)(b)和 P (ε)(c)的電狀態(tài)方程[51],采取原子單位制Fig.3.Electric equations of state of the form ε (D)(a),D (P)(b),and P (ε)(c)in PbTiO3,plotted for fields constrained to lie along the [001],[110],or [111] directions[51].All units are a.u..

考慮原胞形變時,應(yīng)保持材料每個面上的自由電荷不變,即電位移矢量的通量不變,此時可定義約化電位移參數(shù)不 隨原胞形變而變化.可以類似地計算得到開路邊界條件下的應(yīng)力張量.同樣可以用固定電位移方法計算開路邊界條件下的介電張量,其中是電容矩陣的逆,這可以通過有限差分法得到.Jiang 等[52]利用固定電位移方法計算了一些鐵電體的壓電系數(shù).基于開路邊界條件,這一方法也適用于計算模擬鐵電電容器[50]、超晶格體系[53–55]以及金屬-氧化物界面[56–58].基于固定電位移方法,還可以計算饒曲電張量,其中ναβγ=是應(yīng)變梯度張量[59].固定電位移方法目前已在ABINIT 中得到實現(xiàn).

除固定電場法和固定電位移法之外,還有固定電極化法[60].但該方法在物理上對應(yīng)固定極化電荷,實驗上難以實現(xiàn),所以應(yīng)用較少.

3 有限磁場下的DFT 計算

與電場類似,物質(zhì)系統(tǒng)對外部磁場的響應(yīng)通常以磁化強(qiáng)度來描述,磁化強(qiáng)度即單位體積內(nèi)的磁偶極矩大小.在量子理論中,磁化可根據(jù)其來源分為兩種: 軌道磁化和自旋磁化,它們分別源自磁場與電子的軌道自由度和自旋自由度的耦合.因此,通常將與這兩個自由度相耦合的磁場分別稱為“軌道磁場”和“Zeeman 磁場”.值得注意的是,由于晶體場劈裂的作用,軌道磁矩通常很小,所以磁性材料的磁性幾乎全部來自于自旋磁矩;另一方面,軌道自由度對自旋磁矩的貢獻(xiàn)一般來自于自旋軌道耦合(SOC)作用,而這一作用往往很小.基于上述兩個原因,大部分研究在考慮磁場對材料的影響時,通常會忽略軌道磁場而只考慮Zeeman 磁場.

然而,與電場不同的是,無論是Zeeman 磁場還是軌道磁場下的體系依然存在基態(tài),這使得對磁場的處理方式至少在原則上是更嚴(yán)格的.盡管如此,當(dāng)前外場下的第一性原理計算工作仍主要集中于外加電場的情形,對外加磁場的關(guān)注較少.鑒于磁場分別與兩個不同的自由度耦合,本節(jié)接下來將分開討論Zeeman 場和軌道場的影響.除非特別指出,本節(jié)的討論都將忽略SOC 作用.

3.1 Zeeman 磁場和自旋-密度泛函理論

Zeeman 磁場得名于Zeeman 效應(yīng),即原子光譜在外磁場下劈裂的現(xiàn)象,其物理機(jī)制是磁場與原子角動量耦合進(jìn)而打破了原先的能級簡并.而在晶體中磁場主要與自旋角動量耦合,使得自旋向上和自旋向下的電子不再等價,進(jìn)而導(dǎo)致能帶的劈裂.嚴(yán)格考慮自旋的密度泛函理論被稱為自旋-密度泛函理論(SDFT),在20 世紀(jì)70年代由Barth 和Hedin[61,62]提出.由于實際體系都具有自旋,在使用中往往不區(qū)分DFT 和SDFT.

在SDFT 中,電子本征態(tài)需寫成二分量的旋量波函數(shù)的形式Ψnk(r)=(ψnk↑(r),ψnk↓(r))T,其中T表示轉(zhuǎn)置.電子哈密頓量同樣需寫成2×2的 算符矩陣(Vxc)σσ′+μB(σ)σσ′·B,其中μB為玻爾磁子,σ=(σx,σy,σz)為泡利矩陣矢量,下標(biāo)σ為自旋指標(biāo),取值為自旋向上↑或自旋向下↓(默認(rèn)自旋量子化軸為z軸),B為外磁場.顯然Zeeman 項僅與自旋自由度相關(guān)而與空間自由度無關(guān),所以Zeeman 場并不改變體系原有的平移對稱性,這也說明了我們?nèi)园驯菊鲬B(tài)Ψnk寫成布洛赫函數(shù)形式的合理性.此時的能量可以視為 2×2 自旋-密度矩陣的泛函(當(dāng)密度矩陣取基態(tài)密度矩陣時,能量泛函取最小值即基態(tài)能):

所以Zeeman 項可以被自然地吸收進(jìn)勢能項,對(9)式的自洽計算方法與零場時完全一致.

準(zhǔn)確地來說,以上介紹的是非共線版本的SDFT[63–67].共線版本的SDFT 假設(shè)電子自旋僅能朝上或朝下[61,62,68,69],從而哈密頓量和自旋-密度矩陣的非對角元消失,自旋磁矩分布也簡化為m(r)=μB(n↑↑(r)-n↓↓(r)).對于實驗上測得基態(tài)為共線態(tài)的體系,可以直接利用共線SDFT 進(jìn)行計算,從而避免非共線SDFT 帶來的大量計算消耗.歷史上共線SDFT 的發(fā)展也早于非共線SDFT.

施加Zeeman 場的方法當(dāng)前主要用于對多鐵材料的磁電張量的計算.Bousquet 等[70]采用(9)式第一性地研究了磁電張量(此處定義的磁電張量與前文相差一個磁導(dǎo)率μ,即B=μH,且他們的工作中考慮了SOC 和朗德g因子),指出電子貢獻(xiàn)的量級可以與晶格貢獻(xiàn)相比,因而不能被簡單忽略.他們計算了不同大小磁場下 Cr2O3電子極化的變化量(如圖4 所示),從而計算得到體系的磁電張量.在另一個工作中,Bousquet 和Spaldin[71]用同樣的方法計算了施以應(yīng)力下CaMnO3磁電張量電子貢獻(xiàn)的非對角元,發(fā)現(xiàn)其呈現(xiàn)出高度的非線性,這也與群論分析對一階響應(yīng)給出的限制一致.但需要注意的是,通過外加磁場的方式計算磁電張量在數(shù)值上是比較困難的: 由于磁場誘導(dǎo)的原子受力很小,往往需要設(shè)置一個精度很高的收斂條件,如文獻(xiàn)[72] 中設(shè)置的受力收斂條件是 <5 μeV/?.外加Zeeman 場的方法目前在ABINIT 中已得到實現(xiàn)[41].

圖4 Cr2O3 的橫向響應(yīng)貢獻(xiàn)[70] .固定離子響應(yīng) α(el) (空心方塊)的貢獻(xiàn)約為總響應(yīng)的四分之一(實心方塊);響應(yīng)的剩余部分(空心圓)來自于外場下的結(jié)構(gòu)畸變,利用波恩有效電荷計算得到Fig.4.Contributions to the transverse response of Cr2O3 [70].The clamped-ion response,α(el) (open squares)contributes approximately one fourth of the total response (filled circles).The remainder of the response,computed using Born effective charges,is due to structural distortions in the applied field (open circles).

3.2 軌道磁場的密度泛函理論

不同于Zeeman 場直接與自旋耦合,軌道場是以磁矢勢A的形式與正則動量p耦合,此時電子的動能部分為,其中磁矢勢的旋度即為磁場,即?×A=B.這一改變的影響是顯著的: 不同于DFT 可以簡單推廣至SDFT,嚴(yán)格考慮磁矢勢需對DFT 做較大的修改;同時學(xué)界對含磁矢勢的DFT 的關(guān)注和發(fā)展較少,甚至可以說迄今仍處于起步階段.本節(jié)將簡要介紹兩種嚴(yán)格考慮磁矢勢的類密度泛函理論——電流-密度泛函理論(CDFT)[73–77]和磁場-密度泛函理論(BDFT)[78,79].在過去十年里,這兩種理論在小分子層面已經(jīng)得到了一定的發(fā)展和應(yīng)用[80–84].除此之外也有一些含磁矢勢的形式理論,感興趣的讀者可以參看文獻(xiàn)[85].

CDFT 由Vignale 和Rasolt[73]在20 世紀(jì)80年代提出,其最大特點是同時以順磁電流密度jp(r)和電子密度n(r)作為基本變量[86],能量泛函在jp和n為基態(tài)順磁電流密度和基態(tài)密度時取最小值,此時對應(yīng)基態(tài)能.CDFT 的KS 方法同樣要進(jìn)行修改,應(yīng)保證KS 系統(tǒng)給出與真實系統(tǒng)相同的基態(tài)順磁電流密度和基態(tài)密度.KS-CDFT 的能量泛函由下式給出:

其 中AKS(r)=A(r)+Axc(r)和分別是KS 電子感受到的有效磁矢勢和有效勢,和分別是交換關(guān)聯(lián)磁矢勢和交換關(guān)聯(lián)勢,自洽計算時應(yīng)使得n和jp均收斂.注意到KS 系統(tǒng)的有效磁矢勢一般不等于真實系統(tǒng)的磁矢勢,這與前文提到的“密度-極化泛函理論”的情況相似,為了使KS 系統(tǒng)能給出基態(tài)密度以外的物理量,必須引入額外的交換關(guān)聯(lián)勢場.

嚴(yán)格來說,以上介紹的是順磁電流-密度泛函理論,而順磁電流密度并不是規(guī)范不變的,所以歷史上也有人試圖以規(guī)范不變的總電流密度j(r)=jp(r)+n(r)A(r)和電子密度n(r)為基本變量建立“全電流-密度泛函理論”[87,88],但這一理論在數(shù)學(xué)上存在諸多爭議[76,89–91].

BDFT 由Grayce 和Harris[78]在20 世紀(jì)90年代提出,但該理論的后續(xù)發(fā)展極少.與CDFT 不同,BDFT 僅以電子密度n(r)作為基本變量,這樣做的代價是舍棄了通用密度泛函的構(gòu)造,如Grayce和Harris[78]所言,“一個磁場對應(yīng)一個密度泛函理論”.但這個代價對KS 方法的影響可能并不大,我們?nèi)匀豢梢园汛艌龅挠绊憳?biāo)記在交換關(guān)聯(lián)能Exc[n;B]里.KS-BDFT 給出的能量和哈密頓量分別為

以現(xiàn)有水平尚無法判斷CDFT 和BDFT 的優(yōu)劣,但可以肯定的是: 這兩個理論的實際應(yīng)用依賴于合適的磁場下交換關(guān)聯(lián)泛函近似——而這正是目前所缺失的.不過有研究指出meta-GGA 或許具有不錯的前景[82].

3.3 現(xiàn)代軌道磁化理論

在3.2 節(jié)可以看到,軌道磁場下的體系嚴(yán)格具有基態(tài),所以有限軌道場下的第一性計算方法原則上并不需要引入軌道磁化Morb的概念.但考慮到其在核磁共振(NMR)[92]、電子順磁共振(EPR)[93]、量子自旋霍爾效應(yīng)(QSH)[94]、軌道磁電耦合[95–98]等領(lǐng)域具有重要意義,本節(jié)將簡略介紹軌道磁化的理論計算.

軌道磁化的量子力學(xué)定義由21 世紀(jì)建立的“現(xiàn)代軌道磁化理論”給出[99–106]:

其中c為光速,μ為化學(xué)勢,fnk為占據(jù)數(shù).可以發(fā)現(xiàn)該理論與現(xiàn)代電極化理論有許多共通之處,但與電極化不同的是: 軌道磁化是單值的,其依賴于Berry 曲率和哈密頓量;且現(xiàn)代軌道磁化理論適用于有限溫和金屬情形.對于零溫且陳數(shù)為零的系統(tǒng),軌道磁化可以寫為更簡潔的形式[104]:

對于軌道磁化本身(作為零場性質(zhì))的第一性原理計算工作可見文獻(xiàn)[93,107],采用贗勢法計算軌道磁化時應(yīng)考慮贗勢的磁平移對稱性(該對稱性將在下文介紹),這種方法被稱為“含規(guī)范的投影綴加波法 (GIPAW)”[108,109].

不同于前文所述的處理有限電場或Zeeman場的方法,尚未有在能量泛函中直接添加-Morb·B項來做外磁場下第一性原理計算的工作,這可能是基于對軌道磁矩貢獻(xiàn)很小的預(yù)設(shè).然而,必須強(qiáng)調(diào)的是,在某些特定材料中,軌道磁矩的貢獻(xiàn)可能與自旋磁矩相當(dāng),甚至有可能占據(jù)主導(dǎo)地位[110,111].

3.4 含磁矢勢的周期性體系計算

本節(jié)介紹的含磁矢勢的周期性體系計算方法來自于Cai 等[112,113]的工作,據(jù)知,這是目前唯一涉及有限軌道磁場第一性原理計算的工作.引入磁矢勢的主要困難依然在于它破缺了平移對稱性,從而導(dǎo)致布洛赫定理失效,所以不能方便地選用平面波基矢進(jìn)行計算.但幸運的是勻強(qiáng)磁場仍保留有一定的對稱性,能給出所謂的“磁布洛赫定理”,該方法正是利用了此性質(zhì).

為方便起見,本段的討論僅限于n0=1,Γ點和正交晶胞的情形,此時ψnk(r)=unk(r),取a1=,磁矢勢取朗道規(guī)范并略去本征態(tài)的下標(biāo).更一般的討論見文獻(xiàn)[113].在這種選取下,ψ(r)在z方向是周期的,并與x,y無關(guān),可以利用傳統(tǒng)的平面波展開處理,故略去此分量.剩余哈密頓量可拆為三部分:

其中Gy=2π/b,Ky=nyGy是y方向的傅里葉指數(shù),稱 (x,Ky)所在的空間為“介空間”.注意到(18)式事實上給出了螺旋線的拓?fù)?如圖5 所示),因而可以定義“弧長參數(shù)”x:=x+aKy/Gy,將表達(dá)式由二維降至一維,f(x):=f(x,Ky)=f(xa,Ky+Gy).再對f(x)進(jìn)行傅里葉變換即可得到倒空間函數(shù)(注意該倒空間是介空間的“倒”空間).這樣做的好處是和分別在倒空間和介空間對角(如同平面波基矢下動能項在倒空間對角):所以和應(yīng)當(dāng)分別在倒空間和介空間中計算[112].勢能依然由實空間直接相乘得到.

圖5 實空間波函數(shù) ψ (x,y)通過兩次傅里葉變換到倒空間函數(shù) c [112] (a)當(dāng) B=0 時,f(x,Ky)可被視為一系列一維周 期函數(shù)或圓環(huán);(b)當(dāng) B=2π/(ab)時,MPBC 使其變?yōu)橐粭l長螺旋線.由此介空間和倒空間內(nèi)的波函數(shù)可等效為一維函數(shù)Fig.5.The real-space wave function ψ(x,y)can be Fourier transformed into reciprocal space c in two steps[112]: (a)At B=0,f(x,Ky)can be regarded as a set of one-dimensional periodic functions,or rings;(b)at B=2π/(ab),MPBC requires to be a long spiral.The resulting wavefunction in intermediate and reciprocal space is effectively one dimensional.

Cai 等[112,113]利用這種方法計算了磁場下的朗道能級、量子阱中的單電子和雙電子能級、氫原子和氫分子能級以及致密氘流體的電子結(jié)構(gòu).對于致密氘流體,他們發(fā)現(xiàn): 在施加強(qiáng)磁場前后,最高占據(jù)分子軌道(HOMO)和最低未占據(jù)分子軌道(LOMO)變化顯著,然而總電荷密度并未發(fā)生明顯改變,如圖6 所示.

圖6 致密氘流體在零場和強(qiáng)場下的電子結(jié)構(gòu)[112] (a)總電荷密度在B 從0 升到104 T 時基本保持一致;(b)B=0(藍(lán)色)和 B=104 T (紅色)時HOMO 態(tài)在不同原子上的電荷密度分布Fig.6.Electronic structure of dense deuterium fluid under zero and intense magnetic fields[112]: (a)Total charge density remains essentially the same as B goes from 0 to 104 T;(b)the charge densities of the HOMO state for B=0(blue)and B=104 T (red)are distributed on different atoms.

4 處理外場的第一性原理有效哈密頓量方法

當(dāng)代的第一性原理計算軟件已經(jīng)發(fā)展得較為成熟,用戶通常只需要輸入少許參數(shù)即可自動化完成計算并輸出最終結(jié)果,如總能量、極化強(qiáng)度、磁矩等,且結(jié)果通常被認(rèn)為是精確可靠的.但第一性原理在高精度的同時也對算力有較高的要求.受限于計算量,第一性方法往往難以便捷地處理大尺度體系,從而限制了對體系部分動力學(xué)和熱力學(xué)性質(zhì)的研究;此外,這種自動化計算的流程近乎于“黑箱”,在未經(jīng)進(jìn)一步分析的情況下其計算結(jié)果并不能帶來明確的物理圖像.有效哈密頓量方法是為了平衡精度和計算量而誕生的,其旨在使用有限數(shù)量的主要自由度來描述原子哈密頓量或其勢能部分(通常稱為“勢能面”),具體參數(shù)則根據(jù)DFT 計算結(jié)果來確定.從統(tǒng)計物理的角度來說,有效哈密頓量方法是通過積掉配分函數(shù)中的高能或無關(guān)的自由度,從而得到了描述低能物理的有效模型.由于其相對較高的精度與相對較低的計算量,有效哈密頓量方法能處理比DFT 更大尺度的體系,從而更便利地研究體系的相變過程、熱力學(xué)性質(zhì)等,因此在磁性、鐵電、多鐵等領(lǐng)域得到了廣泛的應(yīng)用[116].

4.1 磁性體系的有效哈密頓量方法

由于材料磁性幾乎都由自旋磁矩引起,構(gòu)造磁性有效哈密頓量時一般僅考慮自旋自由度,所以也稱之為“有效自旋哈密頓量”模型.同時為簡化討論與計算,我們常將原子的局域磁矩和自旋視為定長的經(jīng)典歐氏矢量而非量子力學(xué)算符,即采用了“剛性自旋旋轉(zhuǎn)近似”,這在自旋較大時是合理的,對該近似的進(jìn)一步討論見引文[117].一旦得到了磁性有效哈密頓量,便可以輕易計算出磁構(gòu)型對應(yīng)的能量,并借助蒙特卡羅(MC)模擬[118]或者利用Landau-Lifshitz-Gilbert(LLG)方程進(jìn)行自旋動力學(xué)模擬[119–122]來確定磁基態(tài),這兩種方法都適用于有限溫度的場景.若更進(jìn)一步地考慮原子位移自由度,有效自旋模型也可以用來做自旋-晶格動力學(xué)模擬[123,124].

有效自旋模型的精度依賴于第一性原理計算的精度,我們將首先介紹磁性體系第一性原理計算中常使用的兩種方法:“DFT+U”方法和約束磁矩方法.由于傳統(tǒng)的交換關(guān)聯(lián)近似很難描述磁性體系中局域性較高的d 電子和f 電子,所以我們常使用所謂“DFT+U”的方法來補(bǔ)償局域電子的強(qiáng)關(guān)聯(lián)效應(yīng),其中U參數(shù)源自Hubbard 模型中的在位能,具體計算中需要根據(jù)經(jīng)驗或?qū)嶒炛颠x取[125–127].此外我們往往需要第一性原理計算能給出不同的磁構(gòu)型,比如我們希望能計算一個體系所有可能的不同鐵磁、反鐵磁態(tài).雖然原則上KS-DFT 可以給出指定對稱性(如指定總自旋量子數(shù))下的最低能態(tài)[62],但在實際中通常是利用約束總磁矩或約束局域磁矩的方法來得到想要的磁構(gòu)型,具體做法是運用拉格朗日乘子法在能量泛函中添加相應(yīng)的懲罰項[65,128,129]或通過施加局域磁場來約束局域磁矩大小或方向[130,131];其中局域磁矩通常被定義為總磁矩分布在原子附近的積分?I是以I原子為球心的截斷球,其事實上假設(shè)了原子磁矩總局域在原子附近,這對于絕大多數(shù)體系來說都是良好的近似.

無外磁場時多體哈密頓量具有時間反演對稱性(無論是否考慮SOC),有效自旋哈密頓量Hspin也應(yīng)具有這個性質(zhì),即要求有效哈密頓量在自旋全體反向時保持不變,所以模型中應(yīng)當(dāng)僅含有自旋的偶數(shù)次項.對于大多數(shù)體系,我們只需要考慮二階相互作用項即可,其表達(dá)式為[116,132]

其中i,j是原子指標(biāo),Si對應(yīng)于第一性原理計算給出的原子局域磁矩或自旋(假設(shè)原子磁矩總與原子自旋共線),Ai和Jij均為 3×3 的實矩陣,分別稱為單離子各向異性(SIA)矩陣和J 矩陣,它們通常依賴于原子的相對位置.不失一般性地,我們可把SIA矩陣Ai選為對稱無跡陣.J 矩陣可分為三部分,[133,134],其中I為3×3的單位矩陣,為各向同性的海森伯交換相互作用參數(shù),為對稱無跡的各向異性交換相互作用(包括Kitaev相互作用)矩陣[135],為反對稱的Dzyaloshinskii-Moriya (DM)相互作用矩陣[136–138].DM項也常寫成Dij·(Si×Sj)的形式,其中是DM相互作用矢量.不同相互作用傾向于使自旋Si指向不同的方向,體系的磁基態(tài)構(gòu)型依賴于各項的互相競爭.通常體系在零場時的總磁矩取向主要取決于SIA 項和Kitaev 項;當(dāng)DM 相互作用較強(qiáng)時,體系可能形成螺旋態(tài)或斯格明子[139–145].若忽略Kitaev 項和DM 項,二階有效自旋模型將退化為經(jīng)典海森伯模型.高階項中得到較多關(guān)注的是點積形式的四階項和高階手性項,它們在某些材料中也具有顯著效應(yīng)[146,147].

常見的有效自旋模型參數(shù)計算方法大致分為兩大類: 能量映射分析和格林函數(shù)方法.前一大類中常用的方法有擬合法和四態(tài)法.擬合法[148–151]指將有效自旋模型預(yù)測的能量與第一性原理計算結(jié)果一一比較,根據(jù)最小二乘法確定參數(shù)的方法.擬合法原則上可以計算所有類型的相互作用,且可以通過數(shù)據(jù)分析確定參數(shù)的不確定度;但該方法需要預(yù)先生成大量磁構(gòu)型,計算量大.由Xiang 等[152,153]提出的四態(tài)法則是一種計算量較小的能量映射分析方法.在這種方法中,我們需要假設(shè)體系僅包含二階相互作用,可發(fā)現(xiàn)要確定其中任何一個獨立的參數(shù)(可以是矢量或矩陣的一個獨立分量),都只需計算4 種指定自旋構(gòu)型的總能量即可求解(指定自旋構(gòu)型一般通過約束局域磁矩實現(xiàn)),其余參數(shù)的影響都能恰好消掉[132,152,153].該方法的缺點是受高階相互作用影響較大,但特定的高階項仍可以通過巧妙選取自旋構(gòu)型消除[151].近年來,為確定復(fù)雜體系的有效哈密頓量具體形式及各項參數(shù)大小,Li 等[154]基于變量篩選(variable selection)算法發(fā)展了一種新型的能量映射分析方法——機(jī)器學(xué)習(xí)方法構(gòu)造有效哈密頓量(MLMCH)方法,其可以高效而準(zhǔn)確地在諸多備選相互作用形式中挑選出重要的相互作用項,從而建立簡潔準(zhǔn)確的有效哈密頓量;相較于后文將要介紹的神經(jīng)網(wǎng)絡(luò)方法,MLMCH方法能給出解析的有效哈密頓量形式.該方法的可靠性已在有效自旋模型中得到驗證,并集成在PASP (property analysis and simulation package)當(dāng)中[132,154,155].格林函數(shù)方法同樣應(yīng)用廣泛[156–161],該方法主要利用了磁力定理[162]: 對基態(tài)的微擾(這里是兩個原子自旋的微小旋轉(zhuǎn))等于固定的基態(tài)勢能下粒子(這里是電子)能量變化之和.該方法的優(yōu)勢是僅需3 次第一性原理計算即可求得所有的二體二階項參數(shù)和雙二次項(屬于四階項)參數(shù),且僅需用到晶胞而不需要能量映射分析方法中為減小周期性邊界條件影響而使用的較大超胞[163,164];然而,該方法依賴于基組的選擇,對于非磁性原子局域磁矩較大的體系難以準(zhǔn)確描述,且對于明顯偏離參考構(gòu)型的自旋構(gòu)型能量有可能預(yù)測不準(zhǔn).

若將外磁場B納入有效自旋模型,只需要給哈密頓量添上Zeeman 項:

其中g(shù)為朗德g 因子,用于描述軌道磁矩對總磁矩的修正.這一方法在研究磁場對磁性物質(zhì)的調(diào)控方面十分常見.典型的研究如利用其計算磁化強(qiáng)度和居里溫度隨磁場的變化,在考慮溫度時也可以研究磁熱效應(yīng)[165].此外,磁致相變也通常是研究的重點,如以此法計算發(fā)現(xiàn)稀土鐵石榴石在補(bǔ)償溫度附近急劇的全體自旋翻轉(zhuǎn)[166].一類有代表性的研究是用磁場誘導(dǎo)斯格明子態(tài)[167–173]: 如文獻(xiàn)[168] 發(fā)現(xiàn)隨著磁場的增大,低溫下的Pd/Fe/Ir(111)體系發(fā)生了由螺旋態(tài)到斯格明子再到順磁態(tài)的相變;文獻(xiàn)[173] 則更為細(xì)致地研究了CrGe(Se,Te)3Janus單層在不同溫度和磁場下的復(fù)雜相圖,如圖7 所示.CrGe(Se,Te)3同時具有很強(qiáng)的DM 相互作用、阻挫效應(yīng)和較強(qiáng)的面外各向異性,這些特性都有助于斯格明子的穩(wěn)定,因而其被認(rèn)為是研究斯格明子的理想材料[139,173–175].

圖7 CrGe(Se,Te)3 Janus 單層的磁場-溫度相圖[173].相邊界由熱容、磁化率、局域自旋手性決定.這8 個相描述為破碎迷宮疇、斯格明子與嵌套斯格明子合相(I)、迷宮疇(II)、破碎迷宮疇與斯格明子混合相(III)、孤立斯格明子與嵌套斯格明子混合相(IV)、孤立斯格明子(V)、雜化斯格明子相(VI,部分斯格明子合并,部分斯格明子保持分離)、飽和鐵磁態(tài)(VII)、順磁態(tài)(VIII).如圖所示為相III (B=1.8T,T=4.14K)、相IV (B=1.8T,T=4.14K)、相V (B=2.4T,T=4.14K)和 相VI (B=2.4T,T=13.3K)的 代表性自旋結(jié)構(gòu)Fig.7.Magnetic field versus temperature phase diagram of the studied CrGe(Se,Te)3 Janus monolayer[173].The phase boundaries are determined by heat capacity,magnetic susceptibility,local spin chirality,as well as snapshots.The eight phases depicted are as follows: fragmented labyrinth domain,skyrmion and skyrmionium mixed phase (I),labyrinth domain (II),fragmented labyrinth domain and skyrmion mixed phase (III),isolated skyrmion and skyrmionium mixed phase (IV),isolated skyrmion (V),hybrid skyrmion phase (VI,for which some skyrmions merge together but others remain isolated),saturated ferromagnetic state (VII),and paramagnetic state (VIII).Representative spin textures are shown for phase III (B=2.4T,T=4.14K),phase IV (B=1.8T,T=4.14K),phase V(B=2.4T,T=4.14K),and phase VI (B=2.4T,T=13.3K).

4.2 鐵電體系的有效哈密頓量方法

一般認(rèn)為,鐵電相變在宏觀上遵循朗道的對稱性自發(fā)破缺理論,材料由無自發(fā)極化的順電態(tài)(PE)向有自發(fā)極化的鐵電態(tài)(FE)的轉(zhuǎn)變就是體系從高對稱態(tài)“破缺”到了低對稱態(tài).從具體機(jī)制上來看,鐵電相變又可分為位移相變和“有序-無序”相變,前者指原子發(fā)生位移使晶胞從高對稱結(jié)構(gòu)(通常具有中心反演對稱性)變?yōu)闊o中心反演對稱性的低對稱結(jié)構(gòu)并產(chǎn)生了自發(fā)極化;后者指晶胞本身存在自發(fā)極化,但因晶胞的隨機(jī)排布宏觀上不顯極化,宏觀上的自發(fā)極化來自于晶胞排布由無序轉(zhuǎn)為有序的過程.真實材料的鐵電機(jī)制可能介于二者之間.Cochran[176]指出,晶體的位移相變是由那些不穩(wěn)定的聲子模式驅(qū)動的,這些模式被稱為“軟模式”,其頻率呈虛數(shù)[177].在相變點附近,某些聲子之間的非諧性相互作用過強(qiáng)以至于原子事實上偏離了原先的平衡位置,在越過相變點之后聲子模式才重新穩(wěn)定.這些事實表明描述鐵電相變(至少是位移相變)的機(jī)理僅需要用到軟模聲子的自由度.

鐵電的有效哈密頓量方法由Zhong 等[178,179]最早于1994年建立,他們用這種方法成功解釋了鐵電材料BaTiO3的相變機(jī)理,是鐵電研究史上的里程碑.鐵電有效哈密頓量選取的自由度一般是局域應(yīng)變η和若干個“局域模式”u,其中局域應(yīng)變又可以分為均勻應(yīng)變ηH和非均勻應(yīng)變ηI,前者表示晶胞的整體形變,后者表示晶胞角落的局部位移、對應(yīng)于Γ點附近的聲學(xué)聲子;局域模式指的是原胞內(nèi)一些原子的集體運動,不同原胞內(nèi)局域模式的組合可以得到非局域的軟模,這種選取方式被稱為“局域模式近似”[180–182].局域模式自由度考慮到了原子的有限位移,所以原則上可以同時描述位移相變和“無序-有序”相變.更嚴(yán)格的局域模式定義需用晶格Wannier 函數(shù)給出[183],但這一形式應(yīng)用較少,常見的做法仍是將局域模式寫成原子位移疊加的形式.

常見的零場下鐵電有效哈密頓量EFE由五項組成[179]:

分別為局域模式自能、長程偶極-偶極相互作用能、短程軟模相互作用能、彈性勢能和局域模式與應(yīng)變的相互作用能.對于鐵電相變,局域模式往往選為Γ點的光學(xué)聲子模式.有些用于研究合金的鐵電有效模型也考慮了表征同位點原子種類差異的構(gòu)型自由度σ(σ=±1)[184]和原子體積差異的局域應(yīng)變自由度ηloc[185]對能量的影響.此外,由于實際體系中也存在由其他軟模聲子驅(qū)動的反鐵畸變(AFD)相變和反鐵電(AFE)相變,因此許多有效模型也包含了AFD 自由度和AFE 自由度[186–190],在這里仍將其統(tǒng)稱為“鐵電有效哈密頓量”模型.

在鐵電有效哈密頓量中考慮外電場的效應(yīng)一般通過引入局域模式的玻恩有效電荷和固定離子的壓電張量e實現(xiàn),即把極化強(qiáng)度寫為.這樣給出的外電場下鐵電有效模型為[191]

該模型與MC 模擬相結(jié)合,可以用來計算電斯格明子態(tài)[192]、反鐵電相變[189]、電熱系數(shù)αECE=(S為熵)[193–196](如圖8 所示)等,考慮撓曲電效應(yīng)時也能計算撓曲電系數(shù)[197].此外,有效哈密頓量中可以比較方便地引入退極化場,這在部分工作中也已經(jīng)得到體現(xiàn)[198,199].另一種常見的方法是利用有效哈密頓量做分子動力學(xué)模擬(MD),目前利用有限電場下MD 方法研究得較多的是鐵電相變、電熱效應(yīng)等[200–202].

圖8 鐵電材料PbSc0.5Ta0.5O3的電熱效應(yīng)[195](a)鐵電材料PbSc0.5Ta0.5O3的極化強(qiáng)度P(,T)關(guān)于沿〈111〉方向施加的電場和溫度T的函數(shù);(b)電熱系數(shù)α 在330 K時隨電場E 的函數(shù)關(guān)系.在研究溫度下時使 α 達(dá)到極大值的電場[(αmax)]和固定溫度時使得r〈11〉達(dá)到最大值的電場[(r〈11〉)] 也標(biāo)記在圖(a)中.χ2也在圖(b)中標(biāo)出以與 α 做對比.r〈11〉 是大致沿 [11],[11] 或 [11] 方向的局域偶極矩的比例Fig.8.Electrocaloric effects of ferroelectric PbSc0.5Ta0.5 O3[195]:(a)Polarization P(,T)of as a function of electric field applied along 〈111〉 direction and temperature T ;(b)electrocaloric coefficient α as a function of electric field at 330 K.The electric field for which α exhibits its maximum[ (αmax)] and the electric field at which r〈11〉 exhibits its maximum [ (r〈11〉)] for the investigated temperatures are shown in panel (a).χ2 is shown in panel (b)to compare it with α.r〈11〉 is defined as the percentage of local dipoles lying near [11],[11],or [11] directions.

Vanderbilt 的鐵電有效模型原先是基于施加周期性邊界條件的超胞設(shè)計的.Fu 和Bellaiche[203]則對此模型進(jìn)行了擴(kuò)展,開發(fā)出適用于開放邊界條件的鐵電有效哈密頓量方法.他們使用這種新方法并結(jié)合外加電場來研究鐵電微納點的鐵電性質(zhì),發(fā)現(xiàn)在電場驅(qū)動下,鐵電微納點可以產(chǎn)生顯著的極化效應(yīng).Prosandeev 等[204]隨后把這種方法拓展到非勻強(qiáng)電場的情形,體系對電場的響應(yīng)由每個局域電偶極單獨貢獻(xiàn)-pi·Ei,其中pi和Ei為i位點處的局域電偶極和電場大小,他們利用這種方法計算研究了橫向非均勻電場對電環(huán)矩(electric toroidal moment)的調(diào)控.

4.3 多鐵體系的有效哈密頓量方法

多鐵材料可分為兩類: I 型多鐵和II 型多鐵.I 型多鐵通常是良好的鐵電體,并且鐵磁和鐵電相變的溫度可以遠(yuǎn)高于室溫,但此類材料內(nèi)部的磁性與鐵電性的耦合通常很弱.II 型多鐵是新型的多鐵性物質(zhì),其鐵電性僅存在于磁有序狀態(tài)中,并且是由特定的磁序引起,因此該種材料的磁電耦合較強(qiáng).目前主要有兩種基于有效哈密頓量研究多鐵體系的方法.

其一是基于有效自旋哈密頓量方法.但因為有效自旋模型不顯含電場和極化相關(guān)的信息,該方法往往需要額外的DFT 運算.Sasani 等[205]借助海森伯模型研究了GdFeO3的非線性磁致極化現(xiàn)象并給出了與實驗符合的結(jié)果,他們通過旋轉(zhuǎn)Gd 子晶格的磁序?qū)O化強(qiáng)度的變化量表示為G 型反鐵磁序(是Gd 原子和Fe 原子的主要磁基態(tài))強(qiáng)度的函數(shù),這一過程中需要固定Gd 子晶格的磁序但弛豫原子位置和Fe 子晶格的磁矩,并且每一步都需要利用MTP 計算極化強(qiáng)度.Xu 等[37]將有效自旋模型與固定電場方法結(jié)合,提出了以極化強(qiáng)度和DM相互作用為中介、用外加電場可控可逆地調(diào)節(jié)I 型多鐵體系的磁性拓?fù)浜傻臋C(jī)制(EPDQ 機(jī)制),其中有效自旋模型參數(shù)需根據(jù)電場誘導(dǎo)的體系結(jié)構(gòu)確定.

其二是將有效自旋哈密頓量納入到鐵電模型當(dāng)中來構(gòu)建統(tǒng)一的多鐵有效模型,即:

此時的有效自旋模型Hspin還包含了自旋與其他自由度的耦合項[206–211]以及長程偶極相互作用項[206].

考慮多鐵材料在外電場下的響應(yīng),則此時的極化強(qiáng)度P原則上也應(yīng)當(dāng)是 ({S},{u},{η})的函數(shù):

然而,對于I 型多鐵材料,由于磁性和鐵電性互相較為獨立,可以忽略自旋對極化的貢獻(xiàn),極化強(qiáng)度仍寫為.但Bhattacharjee 等[212]通過MC 模擬證實了電場對極化的切換依然能對體系磁序產(chǎn)生重大的影響,這也與EPDQ 機(jī)制一致.而對于磁電耦合顯著的II 型多鐵,我們必須考慮自旋對極化的影響: 例如Xiang 等[213–216]通過建立自旋序誘導(dǎo)的統(tǒng)一極化模型P({S},{u},{η}),成功解釋了II 型多鐵性材料的物理機(jī)制.II 型多鐵性材料的原子位移通常不大,所以可以把極化強(qiáng)度近似寫為電子貢獻(xiàn)和離子-晶格貢獻(xiàn)之和P=Pe({S})+Pion,latt({u},{η}),其中電子貢獻(xiàn)又可寫為,Pi(Si)是在位貢獻(xiàn),Pij(Si,Sj)是位間貢獻(xiàn),〈i,j〉表示最近鄰.電子貢獻(xiàn)與離子貢獻(xiàn)的系數(shù)都能通過四態(tài)法計算得到[215,216].然而,據(jù)知目前尚未有基于統(tǒng)一極化模型建立外電場下的多鐵有效模型的工作.

將外磁場引入多鐵有效模型的方法與自旋有效模型一致,即僅需考慮(20)式中的Zeeman 項即可.Kornev 等[206]采用這種方法計算了BiFeO3塊體的磁電系數(shù),與實驗符合較為良好;Xu等[116,217]基于此模型系統(tǒng)地研究了反自旋-電流相互作用(可視為一種特殊的DM 相互作用)對BiFeO3磁結(jié)構(gòu)的影響(如圖9 所示).特別地,對于II 型多鐵性材料,還可以利用統(tǒng)一極化模型研究和自旋磁矩相關(guān)的磁電效應(yīng).

圖9 塊體BiFeO3在各種C1和C2值下預(yù)測的磁性結(jié)構(gòu)[211],其中C1和C2分別是第一近鄰和第二近鄰的反自旋-電流相互作用系數(shù) (a)計算得到的相圖與C1和C2 的函數(shù)關(guān)系,藍(lán)色十字標(biāo)志和黑色圓點分別代表來自前人選取的C1值(此時C2=0)和 C2 值(此時 C1=0),藍(lán)色三角表示通過密度泛函理論計算得到的結(jié)果,黑線是磁場大小為18 T 時 [ˉ110] 螺旋相向反鐵磁相轉(zhuǎn)變的臨界相;(b)圖示展示了5 種類型的磁螺旋的傳播方向,對于每種類型,紅色、藍(lán)色和黃色分別代表了不同傳播方向的等效磁螺旋Fig.9.Predicted magnetic structures at various C1 and C2 values for bulk BiFeO3 [211],C1 and C2 are coefficients of inverse spin-current interaction for 1st nearest neighbors and 2nd nearest neighbors,respectively: (a)Calculated phase diagram as functions of C1 and C2,the blue cross symbols and black circles are C1 (with C2=0)or C2 (with C1=0)values from previous studies,respectively,the blue triangle symbols are calculated by density functional theory,the black lines are determined by considering the critical magnetic field of 18 T changing the [ˉ110] cycloid to antiferromagnetism;(b)illustration of the propagation directions of the five types of cycloids,for each type,equivalent cycloids of different propagation directions are shown in red,blue,and yellow colors.

4.4 機(jī)器學(xué)習(xí)勢函數(shù)方法

近十年來,機(jī)器學(xué)習(xí)(ML)方法和材料建模與計算的融合成為發(fā)展最快且最具前景的研究領(lǐng)域.相較于前文討論的傳統(tǒng)方法,機(jī)器學(xué)習(xí)方法在保持高精度的同時顯著提升了計算效率,使得原先需要數(shù)小時甚至數(shù)天的計算在幾秒或更短的時間內(nèi)完成.大部分的材料機(jī)器學(xué)習(xí)模型是基于核的機(jī)器學(xué)習(xí)算法,該算法利用材料的描述符作為輸入,基于核脊回歸(KRR)和高斯過程回歸(GPR)等方法學(xué)習(xí)輸入描述符與相應(yīng)材料能量之間的映射函數(shù)[218–223],如原子間勢場模型DPMD[223].然而,近年來這些方法逐漸被性能更優(yōu)秀的圖神經(jīng)網(wǎng)絡(luò)(GNN)算法取代[224–228].圖神經(jīng)網(wǎng)絡(luò)采用連通圖來表示材料的幾何結(jié)構(gòu),網(wǎng)絡(luò)的圖學(xué)習(xí)結(jié)構(gòu)表示可以直接且自然地從輸入結(jié)構(gòu)中學(xué)習(xí)而無需手動構(gòu)建,因此可以被視為一種端到端學(xué)習(xí)的自然描述符.本節(jié)介紹的機(jī)器學(xué)習(xí)方法特指利用神經(jīng)網(wǎng)絡(luò)、尤其是圖神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)并構(gòu)造有效哈密頓量的方法.

據(jù)知,最早基于GNN 系統(tǒng)性地構(gòu)造磁性有效哈密頓量的工作為Yu 等[229]提出的SpinGNN,該網(wǎng)絡(luò)全面地考慮了原子位移自由度和磁矩自由度,對于共線和非共線磁矩均能適用.SpinGNN由兩套獨立的子網(wǎng)絡(luò)組成: 海森伯邊圖神經(jīng)網(wǎng)絡(luò)(HEGNN)和自旋-距離邊圖神經(jīng)網(wǎng)絡(luò)(SEGNN),如圖10 所示.前者使用GNN 中更新的邊特征來映射海森伯相互作用系數(shù),支持學(xué)習(xí)基本的海森伯模型;后者使用原子距離和原子間非共線自旋的點積來初始化網(wǎng)絡(luò)中的邊,可以學(xué)習(xí)高階的磁相互作用;兩套子網(wǎng)絡(luò)可以獨立或共同運行,對于一般的海森伯型,相互作用占絕對主導(dǎo)的材料,只需使用HEGNN,而對于相互作用形式更為復(fù)雜的材料,應(yīng)當(dāng)同時使用兩套網(wǎng)絡(luò),此時體系能量由二者共同給出,Etot=EHEGNN+ESEGNN.經(jīng)測試,該網(wǎng)絡(luò)在預(yù)測多種材料的磁相互作用能和自旋-晶格動力學(xué)模擬方面展現(xiàn)出了優(yōu)越的性能.

圖10 SpinGNN 框架[229],SpinGNN 包含海森伯邊圖神經(jīng)網(wǎng)絡(luò)(HEGNN)和自旋-距離邊圖神經(jīng)網(wǎng)絡(luò)(SEGNN)(a)HEGNN 利用更新后的GNN 邊特征作為Heisenberg 系數(shù),構(gòu)建Heisenberg 型的磁勢;(b)SEGNN 利用自旋-距離邊晶體圖,以自旋矢量的點乘和鍵長初始化邊,可以構(gòu)建一般的高階磁勢,∥ 表示拼接Fig.10.The SpinGNN framework [229],illustration of the SpinGNN including the Heisenberg Edge GNN (HEGNN)and Spin-Distance Edge GNN (SEGNN): (a)HEGNN utilizes the updated edge feature of GNN as the Heisenberg coefficients and builds a Heisenberg-based magnetic potential;(b)SEGNN utilizes the spin-distance edge crystal graph which initializes the edge with the dot product of the spin vector and bond length and builds a high-order general magnetic potential,∥ denotes concatenation.

最新的圖神經(jīng)網(wǎng)絡(luò)技術(shù)是E(3)-等變圖神經(jīng)網(wǎng)絡(luò)(ENN)[230–234],它在原有的圖神經(jīng)網(wǎng)絡(luò)基礎(chǔ)上將體系在三維歐氏群E(3)(平移、旋轉(zhuǎn)、反演)作用下的等變性內(nèi)建在網(wǎng)絡(luò)之中,從而進(jìn)一步減少了訓(xùn)練量并提高了預(yù)測精度.形式上,一個函數(shù)f:X→Y被稱為“G-等變的”指群G的作用與f對易,即:DY[g]f(x)=f(DX[g]x),?g∈G,?x∈X,其中DX[·]和DY[·] 是G在線性空間X和Y上的群表示.在我們考慮的情況下:G是三維歐氏群E(3),X,Y分別表示網(wǎng)絡(luò)的輸入和輸出信息,f是E(3)-等變神經(jīng)網(wǎng)絡(luò)(ENN)實現(xiàn)的映射,其內(nèi)部通常根據(jù)O(3)群的不可約表示對原子特征進(jìn)行直和分解,如NeqiuIP[231]和Allegro[232]這兩個最近提出的E(3)等變的原子間勢場模型.

基于E(3)-ENN 構(gòu)建磁性有效哈密頓量的工作同樣來自于Yu 等[235],他們成功地將E(3)-ENN擴(kuò)展至包含時間反演變換的時間反演-E(3)等變神經(jīng)網(wǎng)絡(luò)(TENN-e3).這一擴(kuò)展是非平凡的,因為原E(3)等變性只考慮到了實空間的對稱變換而沒有涉及自旋空間和速度空間等,具體差異體現(xiàn)在自旋和速度矢量應(yīng)在時間反演下變號而實空間里的矢量、張量在時間反演下不變.TENN-e3 的輸入不僅包含原子間的相對位移rij和原子電荷數(shù)Zi還 包含原子自旋Si,內(nèi)部根據(jù)O(3)群的不可約表示和時間反演下的奇偶性對原子特征進(jìn)行直和分解,因而該方法能適用于考慮或不考慮SOC、共線或非共線的磁性系統(tǒng),其最終輸出任意滿足對稱性要求的標(biāo)量、矢量和張量,比如能量和KS-SDFT 哈密頓量,因此也能自然推廣至外加Zeeman 場的情形.

當(dāng)前使用圖神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)外電場下的有效哈密頓量的研究仍處于初級階段,主流方法仍是基于描述符.Ma 等[236]利用DPMD 建立了以往模型難以描述的鐵電材料HfO2的勢場,并采用外電場下的MD 方法模擬了鐵電反轉(zhuǎn)驅(qū)動的氧離子輸運.但他們的模型并不能預(yù)測玻恩有效電荷和極化強(qiáng)度,模型在外電場下的響應(yīng)仍需通過傳統(tǒng)方法計算.Zhang 等[237]利用深度偶極子模型成功預(yù)測了原子偶極矩以及絕緣體的介電響應(yīng);然而在這種方法中,玻恩有效電荷是通過預(yù)測偶極矩間接計算得到的,可能存在精度偏差.此外,該方法需要同時運用兩種模型進(jìn)行電場作用下的MD 模擬,因此所需計算成本約為原方法的2 倍.

5 總結(jié)與展望

本文的目的是介紹有限外電場和外磁場下的DFT 計算方法和有效哈密頓量模型方法.首先,回顧了現(xiàn)代電極化理論,以及兩種基于此理論構(gòu)建的有限電場下的DFT 算法.然后分別討論了含Zeeman磁場和軌道磁場的密度泛函理論以及與之相關(guān)的現(xiàn)有計算方法.隨后,我們的關(guān)注點轉(zhuǎn)向了鐵電體系和磁性體系的傳統(tǒng)有效哈密頓量方法,以及這些方法在外場下的擴(kuò)展和應(yīng)用.最后介紹了神經(jīng)網(wǎng)絡(luò)在處理外場下周期性體系的應(yīng)用和發(fā)展情況,我們堅信這項技術(shù)在未來將會有更多突破和發(fā)展.

通過對當(dāng)前研究的全面回顧,我們發(fā)現(xiàn),當(dāng)前的外場下DFT 計算方法仍存在諸多挑戰(zhàn)和不足.有限電場方法受制于k點采樣的限制,且無法突破布洛赫定理的框架;有限軌道磁場的計算方法尚處于起步階段,主要由于缺乏適宜的基函數(shù)和交換關(guān)聯(lián)泛函近似.我們預(yù)見,外場下DFT 方法的后續(xù)研究可能會集中在開發(fā)新的基函數(shù)和交換關(guān)聯(lián)泛函近似,同時贗勢對于外場響應(yīng)的研究也可能成為未來的研究焦點.值得我們注意的是,靜態(tài)DFT并非是研究外場對物質(zhì)作用的唯一方法.含時密度泛函理論(TDDFT)[238,239],尤其是含時電流-密度泛函理論(TDCDFT)[240],可能提供了一個更為自然的框架來描述外場與物質(zhì)的相互作用.過去十年中,TDDFT 已經(jīng)取得了顯著的進(jìn)展,我們預(yù)期TDDFT和TDCDFT 的進(jìn)一步發(fā)展及其與機(jī)器學(xué)習(xí)技術(shù)的結(jié)合將是研究周期性體系在外加電磁場下行為的有效工具.