面向鐵電相變的機器學習: 基于圖卷積神經網絡的分子動力學模擬*

歐陽鑫健 張巖星 王之龍 張鋒 陳韋嘉 莊園 揭曉 劉來君 王大威?

1) (西安交通大學電子與信息學部微電子學院,西安 710049)

2) (河南師范大學物理學院,新鄉 453007)

3) (桂林理工大學材料科學與工程學院,桂林 541004)

鐵電材料廣泛應用于功能器件中,對鐵電體進行方便、準確的理論建模,是一個長期被關注的問題.本文提出了一種基于圖卷積神經網絡的鐵電相變模擬方法,利用圖卷積神經網絡對鐵電材料的勢能面進行原子層面的建模,并將得到的神經網絡勢函數作為計算器,以驅動大體系的分子動力學模擬.給定原子位置,訓練好的圖卷積神經網絡能夠給出勢能的高精度預測,達到每原子1 meV 級別,與從頭算(ab inito)精度基本相當,同時在計算速度上相比從頭算方法有數個數量級的提升.得益于神經網絡的高精度和快速預測能力,結合分子動力學模擬,本文對兩種不同類型的鐵電材料——GeTe 和CsSnI3 進行研究,成功模擬了它們隨溫度發生的結構相變,模擬結果和實驗相符合.這些結果說明了圖卷積神經網絡在鐵電體建模和相變模擬應用中的準確性和可靠性,為鐵電體的理論探索提供了一個通用建模方法.

1 引言

鐵電體是指具有鐵電效應的一類材料,經過一個世紀的發展,鐵電材料在現代電子元器件中得到了廣泛的應用[1-4].在鐵電材料的理論研究和數值模擬中,一個基礎性的問題是對鐵電體的能量和原子受力情況進行精確建模,以便得到一個準確的勢函數和力場模型,從而進行大體系的模擬.傳統的建模方法主要包括從頭算(ab initio)方法和物理建模[5-8].從頭算方法利用密度泛函理論(density functional theory,DFT),從第一性原理出發,通過求解多體薛定諤方程得到體系的能量和力學性質,其計算精度高,但是由于計算成本大、效率低,往往只能用于較小的體系和小步數的模擬.

針對上述問題,一種常用的方法是對一個給定的體系進行物理建模.物理建模通過提取體系的物理特征,利用少量的參數對體系的能量給出公式，進行盡量準確的描述.公式中的參數需要通過從頭算或實驗的結果來進行擬合,因此要求研究者對系統有比較深刻的理解,并且具備一定的調參經驗.常用的物理模型包括有效哈密頓量模型[5,6]、核-殼模型[7,8]等.受限于有限的參數,這些模型的精度往往受到限制,同時對于各種各樣的材料結構不具備可遷移性,在應用于復雜體系時亦需要經過繁瑣的人為調參過程.隨著近年來機器學習勢函數[9-14]的發展,其較高的精度和良好的通用性為解決上述問題提供了一個可能的方案.

不同于傳統的以模型為核心的建模方法(常被稱為模型驅動的方法),機器學習方法是一種以數據為核心(數據驅動)的建模方法,能夠從大量數據中自主學習.機器學習中的一個重要研究領域是以神經網絡為基礎的深度學習方法[15].在深度學習中,通過設計合理的數據集和網絡結構,研究者能夠方便地對復雜的體系進行建模,除了需要設置訓練的超參數外,并不需要較多的人為干預.同時,機器學習方法的精度較高,并且能夠容易地遷移到其他體系中.

在眾多的機器學習方法中,近年發展的圖卷積神經網絡[16-20]具有顯著的特點,能夠很好地用于對材料的勢能面和力場進行建模.圖卷積神經網絡將原子間的幾何關系映射成一張無向圖,圖里的結點表示原子,邊表示具有相互作用的原子對,同時對結點和邊賦予相應的特征.Gilmer 等[16]通過消息傳遞機制重新解釋了圖卷積神經網絡的工作原理,因此圖卷積神經網絡也常常被稱為消息傳遞神經網絡.這里,“消息”指的是原子種類特征和原子位置特征(如原子間相對距離、角度等)在高維空間的抽象表示,通過消息在圖上結點之間的傳遞來學習原子間的相互作用關系,進而預測材料所處的能量狀態和原子受力等性質.作為一種數據驅動的機器學習方法,圖卷積神經網絡具有很強的學習能力并能達到較高的預測精度[19].

值得指出的是,盡管基于機器學習方法適用面很廣,原則上可以對任意材料建立力場,并進行分子動力學模擬.但在實踐當中,機器學習方法的應用還是應當針對具體問題進行優化,以期獲得最好的效果.本研究聚焦鐵電材料的相變模擬,相比于對分子體系和對材料生長等其他類型的模擬,這一類工作具有以下幾個顯著的特點.

首先,本文所使用的圖卷積神經網絡最初是針對分子體系的模擬開發的[21,22].為了對鐵電材料相變進行模擬,需要對其進行修改,將其應用拓展到具有周期性條件的晶體,以適用于對鐵電材料的建模.這也說明圖卷積神經網絡方法雖然具有通用性,但在研究中仍需結合相應的體系對其框架進行調整.

其次,鐵電材料相變過程通常涉及復雜的原子位移及其導致的相結構變化,但是其晶格形變較小,很少出現極大應變的情況(多數情況不會大于5%),會保持準立方的結構.因此,在研究當中不需要考慮外加極大應變的情況,并且準立方結構的保持使得構型空間大大減小.這一類模擬與尋找材料在超高壓情況下新類型的結構相是不同的.這些因素應當在研究當中予以充分考慮才能更好地實現對鐵電相變的模擬.

最后,鐵電相變過程中不同相之間能量差別往往很小,對于機器學習方法的預測精度要求較高(誤差應當達到每原子1 meV 的量級).因此,在相變模擬的研究中,要特別關注所訓練的圖神經網絡的適用范圍,使其盡可能在低能部分獲得高精度的預測效果.在樣本生成中體現為需要特別關心低能樣本,特別是與鐵電雙勢阱所在的能量范圍內的樣本.認真考慮上述因素并對所使用的圖神經網絡進行適當調整,才能夠比較好地實現對鐵電相變的模擬.

本文利用圖卷積神經網絡對鐵電材料的勢能面和力場進行建模,并結合分子動力學模擬,研究其隨溫度變化的性質.具體說來,以單硫族化合物GeTe 塊體和鹵族鈣鈦礦CsSnI3兩種具有一定代表性的鐵電材料為研究對象,成功復現了GeTe和CsSnI3隨溫度變化而發生的結構相變,預測了他們的相變序列和相變溫度.值得注意的是,對于GeTe,模擬得到了其在相變附近反常的負熱膨脹效應[23,24];對于CsSnI3,準確獲得了相變過程中伴隨的SnI6八面體轉動模式的變化[25,26].對于這些細微的結構變化,模擬結果能夠和實驗符合,驗證了圖卷積神經網絡在鐵電材料建模中的準確性和可靠性.

2 方 法

為了實現圖卷積神經網絡方法對鐵電相變的模擬,有幾個必要的步驟,總體的流程如圖1 所示,包括數據集的生成、圖卷積神經網絡的搭建以及分子動力學模擬三個部分.

圖1 方法流程圖,包括從頭算分子動力學(ab inito molecular dynamics,AIMD)采樣、圖卷積神經網絡(graph convolutional neural network,GCN)搭建和分子動力學模擬三個部分Fig.1.Workflow of this study,including ab inito molecular dynamics (AIMD) sampling,graph convolutional neural network (GCN) construction and MD simulations.

從頭算分子動力學(ab initomolecular dynamics,AIMD)用于生成訓練數據集,利用分子動力學能夠幫助產生分布合理、符合原子運動規律的樣本.隨后利用生成的樣本集訓練圖卷積神經網絡,可以得到一個精確的勢函數F:{R,Z}→E,F,σ.這一勢函數的輸入是體系的位置信息{R}和原子種類信息(即原子序數Z),輸出是體系的能量E、各個原子的受力F和體系所受到的應力σ.

通過編寫神經網絡勢函數和分子動力學模擬軟件的接口,將圖卷積神經網絡作為分子動力學模擬的計算器,提供每一步模擬所需要的體系性質(即能量、原子受力和體系應力),從而驅動分子動力學的模擬過程.在分子動力學模擬中,通過設置一定的條件(如溫度、壓強等)來模擬體系狀態隨時間發生的變化,生成一系列軌跡.對軌跡進行分析便得到體系在對應條件下的結構性質,如體系的原子位移、自發極化PS和相變溫度TC等.

2.1 圖卷積神經網絡

圖卷積神經網絡代表一類基于圖拓撲結構的網絡框架.在材料建模中,圖卷積神經網絡通常由消息編碼模塊、相互作用模塊和輸出模塊3 大模塊組成.不同的圖卷積神經網絡模型在各模塊的實現方式上有所差別,但其中的關鍵思想是一致的: 通過一組基函數和消息編碼模塊提取原子的種類信息和位置信息,將其映射為更高維的向量(即在圖上傳遞的“消息”),繼而利用相互作用模塊學習原子間的相互作用,最終輸出模塊將消息映射為體系的總能量.

本文使用的圖卷積神經網絡模型為DimeNet++[21,22],如圖2 所示.DimeNet++是近年來圖卷積神經網絡的一個成熟模型,最初是針對分子開發的,本文將其改進并拓展到具有周期性條件的晶體中,以適用于鐵電材料的建模.具體來說,在尋找原子的近鄰原子和原子間距離的計算中,考慮了晶體的周期性;由于鐵電晶體中原子的周圍環境沒有分子體系豐富,相比原來的DimeNet++,使用更少的基函數數目和相互作用層(5 個徑向基函數(radial basis functions,RBF)和3 個角向基函數(sphere basis functions,SBF),2 層相互作用層)已經能夠較好地表征原子的周圍環境.圖卷積神經網絡的各個模塊和文獻[19,22]保持一致,本文著重于闡述圖卷積神經網絡的設計思路和整體框架.

圖2 (a)圖卷積神經網絡框架,以改進后的DimeNet++為例.各個模塊的具體結構和文獻[22]一致;(b)相互作用模塊,包括消息傳遞 finter 和消息更新 fupdate 兩個過程Fig.2.(a)Architecture of the GCN model,a refined Dime-Net++,where the design of blocks are inherited from Reference [22];(b) interaction blocks,including message interaction and message update functions.

相互作用模塊中包含消息傳遞和消息更新兩個過程(圖2(b)),通過消息的傳遞和更新來學習原子間的相互作用.消息傳遞是在圖拓撲結構中的結點之間進行的,每個結點對應一個原子,結點間通過邊連接,邊表示能夠進行消息傳遞的原子對,即能夠發生相互作用的原子對.同時還通過截斷半徑定義了原子之間的作用范圍,從而限制了能和中心原子直接發生相互作用的鄰近原子數目,降低了計算復雜度.

輸出模塊將消息mji進行聚合,并映射為單個原子i的能量Ei.最后一步將所有原子能量相加,得到體系總能量E.消息每經過一次相互作用模塊后,便通過輸出模塊得到第t次相互作用后系統的總能量.模型最終輸出的總能量E是所有輸出模塊輸出的總能量E(t) 的和,即

各種圖卷積神經網絡模型的主要區別在于相互作用模塊,本文使用的相互作用模塊的數學形式如(1)式所示,包括消息傳遞函數finter和消息更新函數fupdate兩個部分.finter和fupdate通過神經網絡實現.消息傳遞函數finter的一個關鍵步驟是針對消息的卷積操作,卷積操作保證了神經網絡的旋轉不變性,同時類比了原子之間的相互作用過程.消息更新函數fupdate借鑒了殘差網絡[28]的思想,利用短路連接來實現消息的更新.相較于圖像識別領域的卷積神經網絡,由于圖卷積神經網絡中圖的拓撲結構往往不是規整的,圖的大小和結點之間的連接關系也隨材料體系的大小和原子的位置發生變化,因此圖卷積中卷積核的大小需要根據近鄰原子的數量發生變化,其卷積的定義亦有所不同(詳見文獻[17]).

需要指出的是,對于某一中心原子j的兩個近鄰原子i和k(即i,k ∈Nj),即使原子i和k之間沒有邊連接,但他們仍然能夠通過和中心原子j的消息傳遞過程發生相互作用.具體來說,在第t次相互作用中,通過對消息的聚合及對的更新來實現消息從原子k到原子i的傳遞.從而每經過一次相互作用模塊,中心原子便可以接收到更遠范圍的消息,和截斷半徑外的原子間接發生相互作用.這是圖卷積神經網絡相較于其他機器學習勢的一個顯著區別.

神經網絡勢函數的一個重要假設是原子相互作用的局域性[9,10],即通過設置合理的截斷半徑來限制原子之間的相互作用范圍,只有在此范圍內的原子才能發生直接相互作用.基于局域性假設,體系的總能量便可以分解為單個原子在其周圍化學環境作用下能量的總和,即,其中Ei是第i個原子處于某一局域環境中所具有的能量.另一方面,局域性假設忽略了原子間長程相互作用的影響,這使得神經網絡勢函數的計算復雜度隨體系的原子數呈線性變化(計算復雜度為O(N)),大大提高了網絡的預測速度.

值得注意的是,雖然局域性假設忽略了原子間長程相互作用,但圖卷積神經網絡中逐次的消息傳遞過程使得原子能夠感知到其截斷半徑以外的原子信息,這與截斷半徑的大小、相互作用模塊的數目有關.截斷半徑越大,相互作用模塊越多,神經網絡能包含的長程相互作用越多,而長程相互作用在大體系的模擬中的作用不可忽視.換言之,更大的截斷半徑和更多的相互作用模塊在一定程度上彌補了神經網絡局域性假設的缺陷,提升了模型的精確度.然而,這也使得圖卷積神經網絡的計算復雜度隨著相互作用次數的增多而有所增加[29].

圖卷積神經網絡是體系能量關于原子位置和原子序數的函數,力和應力則可以通過神經網絡所預測的能量對位置的自動微分[30]計算得到.因此,這里的圖卷積神經網絡應滿足幾個基本的對稱性要求,即能量相對于體系結構的平移不變性、旋轉對稱性和同種類原子的可交換性,這也是所有的機器學習勢應滿足的基本條件.在圖卷積神經網絡中,這些不變性是通過引入相對距離矩陣和圖的卷積操作來實現的.除了對體系勢能面進行建模外,圖神經網絡也能用于預測體系的極化等矢量性質,這就要求網絡相應地具有等變性,即需要使用等變圖卷積神經網絡[27,31].

2.2 數據集的生成

在訓練模型之前,需要通過AIMD 方法生成數據集,即利用從頭算方法,通過分子動力學采樣得到一系列高精度的DFT 數據集.分子動力學模擬是一種常用的模擬方法,能夠從原子或分子尺度研究物質的結構、動力學行為和物理性質[32].它基于牛頓力學和統計力學原理,通過數值計算來追蹤、模擬分子的運動和相互作用,進而得到體系的能量、原子平衡位置、自發極化等性質.

對于GeTe 和CsSnI3,DFT 計算基于綴加平面波方法(projector-augmented wave,PAW)[33,34],在GPAW[35]軟件里實現.樣本體系大小為2×2×2,分別包括64 和40 個原子.計算時使用的交換關聯泛函都是Perdew-Burke-Ernzerhof (PBE)[36]泛函,截斷能分別設為900 eV 和750 eV,k點網格分別為4×4×4和3×3×3 .DFT 計算的收斂標準為GPAW 的默認設置,包括能量、態密度、本征態3 個指標,對應的收斂值分別為5.0×10-4eV/(valence electron),1.0×10-4electron/(valence electron),4.0×10-8eV2/(valence electron),這已經能使CsSnI3和GeTe 的DFT 自洽計算過程中總能量變化收斂至1.0×10-5eV 以內,達到較高的精度.經測試,上述設置對GeTe和CsSnI3的DFT 計算收斂.

分子動力學采樣使用NPT 系綜(具有恒定的粒子數、溫度、壓強),通過Nose-Hoover 算法[37,38]來保持體系的溫度和壓強,外壓強設置為0.分子動力學模擬程序在ASE[39]軟件中實現.對于GeTe,分子動力學模擬的溫度范圍為10—800 K,共產生3202 個樣本;對于CsSnI3,溫度在50—800 K,共產生7908 個樣本.具體而言,在GeTe 和CsSnI3的AIMD 采樣中,為了采集到不同相下的相關樣本,分別通過從高溫的立方相降溫和低溫的基態鐵電相升溫兩個過程來收集樣本.注意在采樣時分子動力學模擬的時間步長應適當大,以減少相似或冗余的樣本構型,本文設置的步長為5 fs.

在本文的分子動力學采樣中,溫度的選取和每個溫度的計算步數并沒有一定的規律.采樣的溫度間隔一般為50 K,而在相變溫度附近則會更密集一些,例如25 K,每個溫度的計算步數也不完全相同,大致在100 步到300 步之間.總之,數據集的構造并沒有一定的規律,往往需要根據訓練出的網絡的預測效果來不斷地迭代、添加新的數據點,因此需要一定的經驗.主動學習方法[40]可以解決此問題,該方法能夠更加智能地構造數據集,在逐步迭代中自動拓展數據集包含的相空間范圍,這也是我們未來的工作.本文主要聚焦于圖卷積神經網絡在鐵電相變模擬中的可靠性.

DFT 樣本的能量分布如圖3 所示,可以看到,樣本主要集中在低能區域,這是由于在關于鐵電相變的勢能面構造中,低能樣本比高能樣本更加重要,更多的低能樣本能有助于神經網絡學習到鐵電體不同相之間細微的差別,而高能樣本則是為了確保在分子動力學模擬中,體系在高溫下能保持動力學穩定性,避免出現崩潰的現象.

圖3 (a)GeTe 和(b) CsSnI3 的DFT 數據集的能量分布柱狀圖,以優化后的立方相的能量為能量零點Fig.3.The energy distribution for the data sets of (a)GeTe and (b) CsSnI3 relative to the energy of corresponding cubic phase.

值得說明的是,盡管分子動力學過程中模擬的誤差會隨著時間步積累,但這并不會對樣本的采集和模型的訓練造成很大的影響.一方面,通過AIMD 生成樣本是一個常用的構造數據集的方法,分子動力學模擬用于生成符合牛頓運動定律的、動力學上合理的樣本.另一方面,在數據集的構造中,數據集上每一個樣本點的能量和原子受力是獨立進行計算的,分子動力學過程中逐步積累的誤差并不會包含在網絡的訓練中.因此,分子動力學模擬中的累計誤差雖然會使獲得的樣本分布發生一定的變化,但對后續圖卷積神經網絡的訓練不會造成本質上的困難.

2.3 神經網絡的訓練

利用神經網絡對鐵電體的勢能面進行建模是一種有監督學習方法.為了能夠在網絡訓練的反向傳播過程中更好地優化網絡參數,需要設計一個合理的損失函數.神經網絡要擬合的是能量關于材料中原子構型的函數關系,因此,損失函數應包括預測能量和真實能量(即DFT 計算的理論值)之間的誤差.更進一步,本文在損失函數中還引入物理先驗知識,即通過將力和應力張量的物理公式加入到損失函數中,使神經網絡在訓練過程中能夠更好地學習體系的勢能面和力場.損失函數的形式為

(2)式中,ρ 是一個常數,起到了權衡能量和力、應力張量的預測精度的作用.對于分子動力學模擬來說,力和應力是關鍵性質,因此將ρ 設置為0.01,以提升力和應力張量在網絡訓練中的優先級.能量可由神經網絡直接預測得到,而原子所受力和體系應力張量則分別通過神經網絡對能量的微分[30]得到,如(3)式和(4)式所示.其中,(4)式根據維里定理(Virial theorem)[41,42]推導得到,式中的α 和β 表示笛卡爾方向,V表示體系的體積,表示原子j作用在原子i上的力.

除了損失函數外,訓練過程中的優化器為AdamW[43],學習率設置為0.001,使用分批訓練的方法,每批包含128 個樣本,即batch size 為128,訓練的總步數為3000 步.在訓練保存的模型中,挑選在驗證集(圖4(b)和圖4(d))上表現最好的模型為最佳模型,用于后續的鐵電相變模擬.

圖4 圖卷積神經網絡模型(GCN)關于GeTe ((a),(b))和CsSnI3 ((c),(d))的測試集和驗證集勢阱預測效果,以優化后的立方相結構為能量零點Fig.4.The test and validation results of the trained graph convolutional neural network (GCN) models for GeTe ((a),(b)) and CsSnI3 ((c),(d)) respectively,where the energy of the optimized cubic phase is set as the reference energy.

3 鐵電相變模擬結果

在所研究的鐵電體中,GeTe 是一種典型的單硫族鐵電材料,在晶相和非晶相之間能發生快速、可逆的結構轉變,以及與結構相關的阻抗特性轉變,是數據存儲設備的理想材料,被廣泛應用在DVD、藍光光盤以及相變開關的制造[44,45].而CsSnI3是一種常見的鹵族鈣鈦礦材料,常溫下直接帶隙為1.3 eV,具有無毒且易于制備的優點[46,47].鹵族鈣鈦礦材料具有優良的光電特性,廣泛應用于太陽能電池、發光二極管和光電探測器等領域[48-51].這兩類材料具有明顯的應用背景,也是比較典型的可能發生鐵電相變的材料.

3.1 圖卷積神經網絡模型的精度

對于GeTe 塊體,3136 個樣本用于訓練,66 個樣本用于測試;對于CsSnI3,7408 個樣本用于訓練,500 個樣本用于測試.訓練好的圖卷積神經網絡模型在測試集上的精度如表1 所列,能量、力和應力的誤差用根均方誤差(root-mean-square error,RMSE)衡量,即其中,是目標預測值,T是真實值,Ntest為樣本數量.表1 中的測試結果顯示,圖卷積神經網絡的預測達到了較高的精度.

表1 GeTe 和CsSnI3 的圖卷積神經網絡模型在各自測試集上的精度Table 1.Prediction accuracy of the trained GCN models for GeTe and CsSnI3 on their test data sets.

圖4 展示了圖卷積神經網絡模型在測試集的能量擬合效果(圖4(a)和圖4(c))和在驗證集勢阱上的能量預測效果(圖4(b)和圖4(d)).其中,GeTe的驗證集是在固定立方相晶格的情況下,將Te 原子整體沿[001] 方向逐步移動生成的.CsSnI3的驗證集是固定立方相晶格,使SnI6八面體沿a軸和b軸逐步發生反相轉動,并沿c軸發生同相轉動(a-a-b+,使用Glazer 記號[52]表示)的樣本,通過PyTilting 程序包[53]生成.圖卷積神經網絡模型準確地預測了GeTe 塊體隨Te 原子位移而產生的勢阱以及CsSnI3塊體隨SnI6八面體轉動而產生的勢阱,能量預測的平均絕對誤差分別為1.473 meV/atom 和0.797 meV/atom.

在GeTe 和CsSnI3塊體的鐵電相變模擬中,分子動力學模擬使用NPT 系綜,并且都從高溫下的立方相結構開始降溫進行.分子動力學模擬的步長設置為2 fs,體系受到的壓強設置為0,每個溫度下的模擬步數為50000 步(即100 ps).對于模擬得到的各個溫度下的結構,序參量(如晶格常數、八面體旋轉角等)是以偽立方坐標系為參照進行統計的,并取每個溫度的最后12500 步(即最后25 ps)進行統計平均,得到體系在該溫度下穩定后的序參量.值得注意的是,本文的研究中除了使用平均晶格的信息作為序參量外,還著重獲得了原子位移隨溫度的變化.對于鐵電材料而言,在已知原子位移的情況下,很容易通過現代極化理論(即Berry phase 方法)[54,55]計算其極化.

3.2 GeTe 塊體的相變模擬結果

本節以GeTe 塊體為例,說明圖卷積神經網絡能用于單硫族鐵電材料的勢能面建模和相變模擬.圖5(a)和圖5(b)分別展示了分子動力學模擬時,GeTe 的晶格常數和Ge-Te 的平均相對位移隨溫度的變化,模擬體系的大小為10×10×10(8000 個原子).Ge 原子和Te 原子的平均相對位移根據(5)式計算得到,其中dGe和dTe分別表示Ge 原子和Te 原子相對于初始立方結構的位移,N為Ge 原子或Te 原子的數目.

圖5 GeTe 塊體的相變模擬結果(a)晶格常數隨溫度發生的變化,紅色虛線框表示鐵電相變附近負的熱膨脹效應;(b) Ge 原子和Te 原子兩者沿 x,y,z 三個方向的平均相對位移(虛線)及其模長(黑色實線)、自發極化(紅色實線)隨溫度的變化情況;(c) Ge 原子和Te 原子的平均相對位移隨體系大小的收斂性測試Fig.5.Phase transition simulations for bulk GeTe: (a)The temperature-dependence of simulated lattice parameters,where the red area indicates the negative volumetric thermal expansion of GeTe near the phase transition;(b) the average relative displacements between Ge and Te atoms during MD simulations and the spontaneous polarization;(c) the convergence of simulation with respect to the system size.

圖5(a)的模擬結果顯示,GeTe 在高溫下呈立方相,在約550 K 時,發生了從立方相至菱方相的相變.在最低的模擬溫度50 K 處,GeTe 的晶格常數約為6.066 ?,夾角約為91.6° (即88.4°).實驗上,GeTe 存在兩種晶相,在室溫下呈現低對稱性的菱方相結構(α-GeTe,R3m),隨著溫度升高,最終GeTe 在約720 K 時發生鐵電相變,轉變為高對稱性的NaCl 結構(β-GeTe,)[23,56,57].分子動力學模擬獲得的相變序列和實驗一致,相變溫度比實驗值[23]低170 K.

圖5(b)計算了Te 原子相對于Ge 原子整體上沿[111] 方向發生的位移,在50 K 時達到0.321 ?,和實驗值(約0.3?)[57,58]接近.在模擬得到體系各個相下原子的平均位移后,利用現代極化理論計算體系的極化,自發極化通過Abinit 軟件計算得到[59].圖5(b)中紅色實線為GeTe 的自發極化隨溫度發生的變化,可以看出極化和Ge 原子與Te 原子之間的相對位移大小隨溫度降低而逐漸變大.Ge 原子和Te 原子之間整體沿[111] 方向的相對位移表明體系自發極化沿[111] 方向.模擬得到GeTe 低溫下的自發極化為0.59C·m2,和理論計算結果(約0.60C·m2)[60,61]接近.

有意思的是,圖5(a)中紅色虛線框表明,在相變溫度附近區域,出現了晶格常數隨著溫度增大而減小的短暫現象,即體系的體積在相變點附近隨著溫度升高而減小,這和之前的研究結果[23,24,62]一致,表明GeTe 在鐵電相變時會出現負的熱膨脹效應,但在遠離鐵電相變溫度處則表現為正常的熱膨脹效應.Dangi?等[24]的研究表明,GeTe 在相變溫度附近反常的熱膨脹現象和聲子的聲學支與軟橫向光學支之間的耦合有關.模擬結果顯示圖卷積神經網絡能夠捕捉到這一細微的變化,充分說明了圖卷積神經網絡的高精度.

3.3 CsSnI3 塊體的相變模擬結果

對于典型的鹵族鈣鈦礦材料CsSnI3,繼續使用圖卷積神經網絡對其勢能面進行建模,研究其包含八面體轉動的結構相變.模擬體系的大小為12×12×12(8640 個原子).和GeTe 類似,首先對CsSnI3的模擬體系大小進行收斂性測試.測試結果表明,分子動力學模擬在12×12×12 大小的體系下是收斂的.

圖6 展示了分子動力學模擬時,CsSnI3的晶格常數(圖6(a))和八面體轉動模式(圖6(b))隨溫度的變化,從模擬結果可以看出,隨著溫度降低,CsSnI3發生了立方相-四方相-正交相的結構相變,同時伴隨著八面體轉動模式從a0a0a0(立方相,C)依次變為a+b0b0(四方相,T)、a+b-b-(正交相,O),相變溫度分別約為300 K(C-T)和200 K(TO).在立方相下,CsSnI3的八面體轉動為零;在四方相下,CsSnI3八面體沿a軸發生同相轉動,即八面體僅繞一個軸發生同相旋轉,在其他兩個軸上旋轉角度為零;在基態正交相下,CsSnI3八面體具有a+b-b-的轉動模式,即SnI6八面體繞一個軸同相旋轉,在其余兩個軸上反相旋轉,且反相旋轉軸上的旋轉角度相同.

圖6 CsSnI3 塊體的相變模擬結果(a)晶格常數隨溫度發生的變化;(b) CsSnI3 在立方相(C)、四方相(T)和正交相(O)下的八面體轉動情況,OOP (out-of-phase)和IP (inphase)分別表示反相和同相轉動Fig.6.Phase transition simulations for bulk CsSnI3: (a)The temperature-dependence of simulated lattice parameters;(b) the change of SnI6 octahedron tilting pattern during the cubic-tetragonal-orthorhombic (C-T-O) phase transition.

4 討 論

本文提出了一個針對鐵電體勢能面建模的通用框架,采用基于圖卷積神經網絡的分子動力學模擬方法對GeTe 和CsSnI3的鐵電相變進行了研究,本節對模擬結果和這一研究方法進行討論.

4.1 GeTe 塊體

對于單硫族化合物GeTe 塊體隨溫度降低而發生的鐵電相變,在從高溫到低溫的降溫模擬中,體系結構從立方相轉變為菱方相.在菱方相下,GeTe 相較于其立方相(順電相)沿[111] 方向發生形變,同時Te 原子相對Ge 原子沿[111] 方向存在相對位移[58].對于這一過程,在早期曾嘗試通過以公式為核心的方法進行建模[64],但是單硫族化合物中存在特殊的“超鍵 (metavalent bond)”[65],并沒有清晰的軟模概念,因此在利用有效哈密頓量方法對其勢能面建模時存在較大困難[66],而圖卷積神經網絡避免了人為構建能量函數,通過消息傳遞過程自動學習原子之間的相互作用,從而更加方便、高效地對單硫族化合物的勢能面進行建模.從這個意義上說,圖卷積神經網絡可以視為“無需物理先驗知識的有效哈密頓量方法”[67].

張佳樂等[68]通過Behler-Parrinello (BP)[9]神經網絡模型對單硫族的塊體材料以及二維單原子層GeTe 和SnTe 的原子間勢函數進行了建模,并結合蒙特卡羅模擬研究了其相變過程.這一工作并未考慮材料的整體應變,因此預測的相變溫度不能和實驗值進行直接比較.這一類的研究顯示,二維單硫族化合物在單胞厚度下依然具有鐵電性,并且相變溫度較高,具有潛在的應用價值[69].利用本文提出的基于圖卷積神經網絡的分子動力學方法,可以更方便地對二維單硫族鐵電材料的相變過程進行全面模擬.

4.2 CsSnI3 塊體

相比于單硫族化合物GeTe 的單個順電-鐵電相變,鹵族鈣鈦礦的相變過程更加豐富,且伴隨著八面體轉動等有趣而復雜的現象.

在鹵族鈣鈦礦中,一般有三種相變的機制,包括陽離子的位移、八面體的旋轉以及八面體的畸變.其中,八面體的旋轉相變是較為常見的相變類型[70].對于ABX3構型的鈣鈦礦結構,其八面體的旋轉相變大多源于空間效應,一般當A位陽離子的半徑相對較小時八面體有發生旋轉的傾向[26],A位陽離子的電子結構也會影響八面體的旋轉[71].在有機無機雜化的鹵族鈣鈦礦中,有機離子與陰離子框架之間的氫鍵是八面體旋轉的主要原因之一[72].

與常規的鈣鈦礦材料(如BaTiO3,PbTiO3)不同,鹵族鈣鈦礦中的八面體轉動一直是一個難以處理的問題,因為相鄰八面體之間互相鉸連,一個八面體的轉動必然會引起周圍八面體的轉動,而且八面體轉動造成的能量變化是非常微妙的.這些因素使得鹵族鈣鈦礦的理論建模更加復雜.在以往的工作中,Chen 等[73]利用有效哈密頓量法對CsPbI3的相變和碘八面體轉動進行模擬,Jinnouchi 等[74]利用高斯核回歸[75]方法對有機鹵族鈣鈦礦Methylammonium(MA)PbI3的八面體轉動相變進行了模擬.

由于圖卷積神經網絡的高精度,能精確捕捉這類晶體中原子微小位移造成的能量變化,并最終擬合出鹵族鈣鈦礦材料中八面體轉動的勢能面.將擬合好的勢函數應用于分子動力學模擬,成功得到了CsSnI3隨溫度變化發生的結構相變,準確得到了各個相下的八面體轉動模式,顯示了圖卷積神經網絡在鹵族鈣鈦礦勢能面建模中的可靠性.除了預測相變的靜態性質外,這一方法還可用于進一步研究鹵族鈣鈦礦相變過程的動態性質.通過觀察原子位置在模擬過程隨時間發生的變化,能更細致地分析鐵電材料在結構相變中的動態性質,如相變的有序-無序性和離子遷移等,正如我們之前的工作[18]所展示的.

4.3 數據驅動的勢能面建模

對于鐵電材料建模的一類常用方法是有效哈密頓量方法,自有效哈密頓量提出的近三十年來,已成功應用于許多鈣鈦礦鐵電材料的靜態和動態性質的研究中[76,77].這一方法將體系的能量描述為關于軟模的函數,分解為若干個能量項(如偶極子自能、偶極子之間的長(短)程相互作用能等)之和,并結合DFT 計算的結果對各個能量項的參數進行擬合.然而,由于有效哈密頓量方法通過軟模的概念對體系自由度進行了簡化,在精度上存在一定的限制.從對系統進行建模的角度來講,GeTe和CsSnI3這兩類材料都給有效哈密頓量等傳統方法帶來一定的挑戰.

首先,對于GeTe 來講,有效哈密頓量方法是對其進行建模的一個可能辦法.這一方法利用軟模的物理概念提取鐵電體勢能面的關鍵特征,然而,由于單硫族化合物獨特的電子特性,存在介于共價鍵和離子鍵(或金屬鍵)之間的metavalent 鍵[18,78],使用有效哈密頓量方法對其勢能面建模存在困難.具體而言,這類材料與典型鈣鈦礦鐵電體(例如BaTiO3)不同,其中原子間具有很強的非線性相互作用,并且這一非線性相互作用不能通過偶極-偶極間相互作用分離出來,所以難以提取出軟模并用有效哈密頓量方法處理[64,66].因此,目前對于單硫族化合物的勢能面建模缺少一個統一、優美的方法.已有的機器學習工作[68]使用了物理描述符(描述符表征了原子的周圍環境和原子之間的相互作用)加全連接網絡的方式,沒有使用比較先進的神經網絡架構.

其次,對于鹵族鈣鈦礦而言,八面體的轉動會導致細微的能量變化,其勢函數更加復雜,難以通過傳統的方法進行精確建模[18,78].而且,八面體旋轉相變對鹵族鈣鈦礦材料的帶隙有著較大的影響,而帶隙直接影響了半導體材料的光電性能,因此鹵族鈣鈦礦材料中的八面體旋轉相變更值得關注.本文通過圖卷積神經網絡對上述兩種材料的勢能面進行了建模,發現圖卷積神經網絡模型達到了令人滿意的精度,能夠用于分子動力學模擬預測CsSnI3隨溫度發生的結構變化.

對于鐵電相變的模擬,在近幾年間發展出了幾種不同機器學習方法.高斯核回歸方法是另一種常用的機器學習方法,通過對所有訓練樣本的每一個原子都提取出相應的描述符,在預測時將體系中的原子環境表示為已知描述符的線性組合,從而得到體系的能量.隨著訓練樣本的增多,高斯核回歸方法要提取并處理、保存的描述符也越多,模型復雜度和內存消耗也隨之增大,因此其可擴展性可能會出現困難[75].與有效哈密頓量和高斯核回歸等方法相比,圖卷積神經網絡能夠實現原子級別的預測,同時具有更高的分辨率和更高的精度(可達到1 meV/atom),而且其模型復雜度和訓練樣本數沒有必然的聯系,因此能夠很好地用于GeTe 和CsSnI3等復雜鐵電體的建模與模擬當中.

除了圖卷積神經網絡外,以神經網絡為核心的機器學習勢函數還包括深勢模型(deep potential,DP)[79]、Behler-Parrinello (BP)模型等.圖卷積神經網絡與其他網絡的重要區別在于其消息傳遞過程.不斷的消息傳遞過程使得原子最終能感知到遠處原子的影響,從而使原子的作用范圍變為半局域,提升了網絡的表達能力.此外,相較于BP模型,圖卷積神經網絡無需預先人為地構造描述符,而是將原子之間的關系抽象成圖,直接從圖網絡中學習到原子之間的相互作用關系,這使得圖卷積神經網絡能夠方便地遷移到其他鐵電材料的建模中,對于各式各樣的鐵電材料更加通用.

本文通過和結構有關的序參量來表征體系的結構變化,進而確定體系的相變溫度.相變溫度亦可通過熱力學積分的自由能方法[80]得到.熱力學積分方法將體系自由能表示為耦合參數λ 的函數,通過對兩個狀態間關于λ 的路徑積分來計算體系的自由能變化[81],而圖卷積神經網絡直接將能量表示為原子種類和坐標的函數.因此,對于基于圖卷積神經網絡的分子動力學模擬來說,直接通過分析體系結構的變化(如晶格常數、原子位移)來判斷相變更為方便.以往的工作[18,19]表明,通過分析體系結構的變化來判斷相變是可靠的.

值得指出的是,本文預測的相變溫度與實驗值有所差異,這和用于DFT 樣本計算時交換關聯泛函的選取有關[82,83].因為圖卷積神經網絡擬合的實際上是DFT 計算的理論值,而DFT 計算的準確度極大地依賴于交換關聯泛函的準確性.以PbTiO3為例,在之前的研究中已經發現,選取不同的交換關聯泛函進行DFT 計算,PbTiO3的相變勢阱,特別是勢阱深度差異較大,導致預測的相變溫度存在較大差別.具體來說,PBE 泛函會導致PbTiO3的勢阱深度比LDA 深很多[84],預測的相變溫度比LDA 高約200 K[19].在不同的交換關聯泛函下,DFT 結果存在不可忽視的差異會導致機器學習預測的相變溫度出現差別.使用更為準確的交換作用勢(例如HSE06 雜化泛函)能夠獲得更好的結果.但在實踐當中,我們發現這類泛函的計算成本太高,應用起來也比較困難.這一方面的問題值得在未來進行更多的探索.

5 結 論

本文利用圖卷積神經網絡對鐵電體勢能面進行建模,為鐵電相變的模擬提供了一種通用的建模方法.以GeTe 塊體和CsSnI3為研究對象,結合分子動力學模擬,成功地模擬了兩者隨溫度發生的結構相變,展示了這一方法的具體應用,希望能為這一方法在更多實際應用中提供有價值的信息.本文中模擬的結果充分說明了圖卷積神經網絡在精度方面的潛力.例如GeTe 在相變溫度附近負的熱膨脹效應和CsSnI3的八面體轉動相變,表明了圖卷積神經網絡在鐵電相變模擬中的可靠性和準確性.

總之,針對鐵電結構相變,經過良好訓練的圖卷積神經網絡能夠達到DFT 的計算精度,其計算成本低,適用于大體系、長步數的模擬.雖然本文聚焦于結構相變的模擬,但是毫無疑問,圖卷積神經網絡方法能夠應用于更廣泛的材料研究中.例如,在未來的工作中還可將嘗試將此方法應用于鐵電材料的熱導率、聲子譜等性質的計算.