外加電磁場下周期性體系的第一性原理計算方法*

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

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

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

3)(人工微結構科學與技術協同創新中心,南京 210093)

1 引言

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

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

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

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

2 有限電場下的DFT 計算

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

2.1 現代電極化理論

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

圖1 兩種原子構成的系統的不同原胞選擇[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 的高明之處在于: 它并不直接定義電極化強度本身,而是定義極化強度在晶體絕熱演化下初末態之間的變化量[18].事實上實驗學家測量的也正是極化強度的變化量而非極化強度本身,如鐵電材料的自發極化是通過測量電滯回線得到的.極化強度的變化量可以被定義為

其中 i和f 分別指代初、末態,總極化強度P可認為是離子極化強度Pion和電子極化強度Pel的和,離子極化強度定義為

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

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

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

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

其中ai是原胞基矢,

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

其中Wannier 函數

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

2.2 固定電場方法和電熱焓

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

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

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

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

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

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

表1 用固定電場方法計算的一些III-V 半導體介電性質與實驗的比較[42],其中玻恩有效電荷張量在材料對稱性下退化為標量,且 d123 定義為 /2,LDA 和PBE 是計算時使用的交換關聯泛函近似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 半導體介電性質與實驗的比較[42],其中玻恩有效電荷張量在材料對稱性下退化為標量,且 d123 定義為 /2,LDA 和PBE 是計算時使用的交換關聯泛函近似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.

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

圖2 固定電壓法和固定電場法的差異[45].當材料發生應變 η時,如果保持電壓?V不變,電場從E=?V/d 變 化為 E=?V/[(1+η)d] ;如果保持電場 E 不變,電壓從 ?V 變為(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.

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

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

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

2.3 固定電位移方法和內能泛函

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

圖3 施加場的方向約束在[001],[110] 或[111] 方向時,PbTiO3 中形式為 ε (D)(ε 為外電場)(a)、D (P)(b)和 P (ε)(c)的電狀態方程[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..

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

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

3 有限磁場下的DFT 計算

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

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

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

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

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

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

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

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

圖4 Cr2O3 的橫向響應貢獻[70] .固定離子響應 α(el) (空心方塊)的貢獻約為總響應的四分之一(實心方塊);響應的剩余部分(空心圓)來自于外場下的結構畸變,利用波恩有效電荷計算得到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,嚴格考慮磁矢勢需對DFT 做較大的修改;同時學界對含磁矢勢的DFT 的關注和發展較少,甚至可以說迄今仍處于起步階段.本節將簡要介紹兩種嚴格考慮磁矢勢的類密度泛函理論——電流-密度泛函理論(CDFT)[73–77]和磁場-密度泛函理論(BDFT)[78,79].在過去十年里,這兩種理論在小分子層面已經得到了一定的發展和應用[80–84].除此之外也有一些含磁矢勢的形式理論,感興趣的讀者可以參看文獻[85].

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

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

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

BDFT 由Grayce 和Harris[78]在20 世紀90年代提出,但該理論的后續發展極少.與CDFT 不同,BDFT 僅以電子密度n(r)作為基本變量,這樣做的代價是舍棄了通用密度泛函的構造,如Grayce和Harris[78]所言,“一個磁場對應一個密度泛函理論”.但這個代價對KS 方法的影響可能并不大,我們仍然可以把磁場的影響標記在交換關聯能Exc[n;B]里.KS-BDFT 給出的能量和哈密頓量分別為

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

3.3 現代軌道磁化理論

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

軌道磁化的量子力學定義由21 世紀建立的“現代軌道磁化理論”給出[99–106]:

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

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

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

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

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

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

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

圖5 實空間波函數 ψ (x,y)通過兩次傅里葉變換到倒空間函數 c [112] (a)當 B=0 時,f(x,Ky)可被視為一系列一維周 期函數或圓環;(b)當 B=2π/(ab)時,MPBC 使其變為一條長螺旋線.由此介空間和倒空間內的波函數可等效為一維函數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]利用這種方法計算了磁場下的朗道能級、量子阱中的單電子和雙電子能級、氫原子和氫分子能級以及致密氘流體的電子結構.對于致密氘流體,他們發現: 在施加強磁場前后,最高占據分子軌道(HOMO)和最低未占據分子軌道(LOMO)變化顯著,然而總電荷密度并未發生明顯改變,如圖6 所示.

圖6 致密氘流體在零場和強場下的電子結構[112] (a)總電荷密度在B 從0 升到104 T 時基本保持一致;(b)B=0(藍色)和 B=104 T (紅色)時HOMO 態在不同原子上的電荷密度分布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 處理外場的第一性原理有效哈密頓量方法

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

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

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

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

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

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

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

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

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

圖7 CrGe(Se,Te)3 Janus 單層的磁場-溫度相圖[173].相邊界由熱容、磁化率、局域自旋手性決定.這8 個相描述為破碎迷宮疇、斯格明子與嵌套斯格明子合相(I)、迷宮疇(II)、破碎迷宮疇與斯格明子混合相(III)、孤立斯格明子與嵌套斯格明子混合相(IV)、孤立斯格明子(V)、雜化斯格明子相(VI,部分斯格明子合并,部分斯格明子保持分離)、飽和鐵磁態(VII)、順磁態(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)的 代表性自旋結構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 鐵電體系的有效哈密頓量方法

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

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

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

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

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

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

圖8 鐵電材料PbSc0.5Ta0.5O3的電熱效應[195](a)鐵電材料PbSc0.5Ta0.5O3的極化強度P(,T)關于沿〈111〉方向施加的電場和溫度T的函數;(b)電熱系數α 在330 K時隨電場E 的函數關系.在研究溫度下時使 α 達到極大值的電場[(αmax)]和固定溫度時使得r〈11〉達到最大值的電場[(r〈11〉)] 也標記在圖(a)中.χ2也在圖(b)中標出以與 α 做對比.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 的鐵電有效模型原先是基于施加周期性邊界條件的超胞設計的.Fu 和Bellaiche[203]則對此模型進行了擴展,開發出適用于開放邊界條件的鐵電有效哈密頓量方法.他們使用這種新方法并結合外加電場來研究鐵電微納點的鐵電性質,發現在電場驅動下,鐵電微納點可以產生顯著的極化效應.Prosandeev 等[204]隨后把這種方法拓展到非勻強電場的情形,體系對電場的響應由每個局域電偶極單獨貢獻-pi·Ei,其中pi和Ei為i位點處的局域電偶極和電場大小,他們利用這種方法計算研究了橫向非均勻電場對電環矩(electric toroidal moment)的調控.

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

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

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

其二是將有效自旋哈密頓量納入到鐵電模型當中來構建統一的多鐵有效模型,即:

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

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

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

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

圖9 塊體BiFeO3在各種C1和C2值下預測的磁性結構[211],其中C1和C2分別是第一近鄰和第二近鄰的反自旋-電流相互作用系數 (a)計算得到的相圖與C1和C2 的函數關系,藍色十字標志和黑色圓點分別代表來自前人選取的C1值(此時C2=0)和 C2 值(此時 C1=0),藍色三角表示通過密度泛函理論計算得到的結果,黑線是磁場大小為18 T 時 [ˉ110] 螺旋相向反鐵磁相轉變的臨界相;(b)圖示展示了5 種類型的磁螺旋的傳播方向,對于每種類型,紅色、藍色和黃色分別代表了不同傳播方向的等效磁螺旋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 機器學習勢函數方法

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

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

圖10 SpinGNN 框架[229],SpinGNN 包含海森伯邊圖神經網絡(HEGNN)和自旋-距離邊圖神經網絡(SEGNN)(a)HEGNN 利用更新后的GNN 邊特征作為Heisenberg 系數,構建Heisenberg 型的磁勢;(b)SEGNN 利用自旋-距離邊晶體圖,以自旋矢量的點乘和鍵長初始化邊,可以構建一般的高階磁勢,∥ 表示拼接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.

最新的圖神經網絡技術是E(3)-等變圖神經網絡(ENN)[230–234],它在原有的圖神經網絡基礎上將體系在三維歐氏群E(3)(平移、旋轉、反演)作用下的等變性內建在網絡之中,從而進一步減少了訓練量并提高了預測精度.形式上,一個函數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分別表示網絡的輸入和輸出信息,f是E(3)-等變神經網絡(ENN)實現的映射,其內部通常根據O(3)群的不可約表示對原子特征進行直和分解,如NeqiuIP[231]和Allegro[232]這兩個最近提出的E(3)等變的原子間勢場模型.

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

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

5 總結與展望

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

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