尖晶石氧化物能量和結構的第一性原理計算和機器學習

李一航 肖 斌 唐宇超 劉 馥 王小夢 劉 軼

(1.上海大學理學院， 上海 200444; 2.上海大學材料基因組工程研究院， 上海 200444)

高通量第一性原理計算方法[1-2]可用于系統研究大量的材料體系和構型， 但龐大的計算量使其成本高、耗時長.機器學習(machine learning， ML)方法[3-4]可以從計算或實驗材料數據中進行數據挖掘或機器學習， 獲得其中隱含的關聯規律， 并用于對未知材料性能的高效預測.本研究基于高通量第一性原理計算數據構建機器學習模型， 可加快新材料的能量和結構預測，指導新材料的實驗合成.

尖晶石材料因其龐大的種類、靈活可調的組分以及優異的物理化學性質而被廣泛應用于鋰電池正極[5-6]、光催化劑材料[7-8]、鐵氧體磁性材料[9]和半金屬材料[10]等領域.尖晶石氧化物可分為正尖晶石結構、反尖晶石結構和混合尖晶石結構[11-13].本研究中的正尖晶石氧化物化學通式為AB2O4， 空間群為Fdm， 其中A， B 陽離子[14]分別處于氧離子形成的四面體間隙(T 位點)和八面體間隙(O 位點)中.通過多種方式， 結合T/O 位點中心原子TC/OC和T/O 位點環境原子TE/OE特征構建的中心-環境(center-environment， CE)特征模型， 如圖1 所示.本研究用73 種元素(見圖2)分別替換尖晶石氧化物AB2O4結構的T 位點和O 位點原子， 獲得了5 329(732=5 329)個可能的尖晶石氧化物.圖2 中的灰色背景元素表示沒有研究， 色標從藍色到粉色作為對應的單原子基態總能量的函數， 正方形和圓形代表機器學習預測為“好”的T 位點和O 位點的穩定性元素.

圖1 正尖晶石慣用單胞和CE 特征模型Fig.1 Conventional cell of the normal spinel type and the CE feature models

圖2 正尖晶石氧化物AB2O4 的組成元素Fig.2 Constitution elements of the normal spinel oxides AB2O4

在計算設計預測新材料時， 要求新材料必須是熱力學穩定或亞穩定的， 而形成能[15-16]的大小是評估新材料是否熱力學穩定的重要判據之一.為了得到準確的能量， 需要對預測結構進行充分弛豫， 弛豫后的計算結構可作為未來實驗驗證結構的理論參考依據.本研究使用高通量第一性原理計算方法研究了5 329 個尖晶石氧化物AB2O4結構的形成能和晶格常數， 提出了一種“中心-環境”(CE)特征模型， 以構建同時包含成分和結構信息的特征， 作為機器學習算法的輸入變量.基于密度泛函理論(density functional theory， DFT)計算數據， 使用隨機森林算法開發了機器學習模型， 對尖晶石氧化物結構的DFT 計算數據進行了學習訓練， 獲得的機器學習模型可以準確有效地預測尖晶石氧化物的形成能和晶格常數.

1 研究方法

1.1 第一性原理計算

本研究采用基于贗勢投影綴加波的第一性原理計算方法， 通過VASP(Vienna Abinitio Simulation Package)軟件包[17]進行計算.交換關聯泛函選取廣義梯度近似-PBE(generalized gradient approximation proposed by Perdew， Burke and Ernzerhof， GGA-PBE)， 平面波截斷能為520 eV， 結構優化過程中保持能量的收斂精度為1×10-5eV/atom， 力的收斂精度為0.02 eV/°A(1°A=0.1 nm).通過收斂性測試， 布里淵區采用Monkhorst-Pack 的7×7×7 網格.本研究中尖晶石氧化物形成能的計算公式為

式中:Espinel是尖晶石晶胞的基態總能量;ET，EO和EO2分別為單個T 位點原子、單個O 位點原子和氧分子的基態總能量;nT，nO和NO分別為晶胞中T 位點原子、O 位點原子和氧原子的數量;N為單胞中的總原子數.

1.2 CE 特征模型

本研究結合元素/單質的基本性質和尖晶石氧化物的組成與結構信息來構造復合特征.每個構成元素選取了元素和相應單質的55 個基本性質[18]， 如原子半徑、氧化態、電負性等(見表1).

表1 元素和單質的基本性質Table 1 Fundamental features of elements and simple substances

基于特征工程理念， 本研究構造了CE 特征模型， 通過將基本屬性投影到由組分和結構信息組成的基組上來構造特征.定義一個由中心原子和周圍環境原子組成的CE 原子集， 其中環境原子包括距離中心原子第1 近鄰到第N近鄰的原子(N=1～10).如果只有一個中心原子及其對應的環境原子， 則稱之為單中心CE(single CE， sCE)特征模型; 如果包括兩個中心原子及其對應的環境原子， 則稱之為雙中心CE(dual CE， dCE)特征模型.

CE 特征模型的基本思想是分別定義區分中心原子和環境原子， 然后將環境原子的基礎性質以與中心原子距離的倒數作為權重進行線性加和， 從而實現將化學成分和結構信息編碼到機器學習特征中.機器學習訓練集D包括特征Fi(i=1～55 為性質序數)和預測目標T， 具體定義為

特征Fi(i=1～55)包括中心原子特征fC，i和環境原子特征fE，i， 即

中心原子的基本性質pC，i可直接作為該中心原子的特征， 即

環境原子的特征是根據每個環境原子j的基本性質pj，i及其到中心原子距離的倒數， 加權求和而成的復合屬性特征， 即

式中，ωj為求和權重， 設環境原子j與中心原子的距離為rj， 則

本研究以單中心和雙中心、標量和矢量等不同方式將中心原子和環境原子特征組合，構建成多種復合CE 特征模型.dCE 特征模型分別選取尖晶石氧化物結構的非等效T 位點和O 位點作為中心.TC和TE代表T 位點的中心原子及其周圍環境， 如圖1(a)中Ⅰ和Ⅱ所示; OC和OE代表O 位點的中心原子及其周圍環境， 如圖1(a)中Ⅲ和Ⅳ所示.CE-1， CE-3， CE-6 和CE-7 為4 個sCE 特征模型; CE-2， CE-4， CE-5 和CE-8 為4 個dCE 特征模型.在這8 個特征模型中， 分別將T 位點或O 位點定義為中心原子， 并將其第1 近鄰至第10 近鄰的原子分別考慮為環境原子.弛豫后的原型鎂鋁尖晶石MgAl2O4幾何構型信息， 被用作構建本研究的尖晶石氧化物結構CE 特征模型時的基礎結構.使用固定的原型結構信息來構造CE 特征氧化物， 可使機器學習模型更易于預測和遷移， 不受個別結構弛豫后幾何構型變化的影響， 從而避免對新未知結構進行昂貴的第一性原理計算.

1.3 機器學習預測形成能和晶格常數

本研究采用隨機森林算法訓練機器學習模型.采用前述的CE 特征模型作為機器學習模型的輸入變量， 預測性能包括尖晶石氧化物結構的形成能和晶格常數等.選取所有數據的80%作為訓練集， 采用5 折交叉驗證， 通過網格搜索確定最佳超參數.剩余的20%數據作為完全獨立的性能測試集， 并不參與模型訓練.將決定系數(coefficient of determination)R2、平均絕對誤差(mean absolute error， MAE)和均方根誤差(root mean squared error， RMSE)作為機器學習模型性能的評價指標， 具體定義分別如下:

式中:yi是實驗或計算獲得的目標值;是單個樣本yi的預測值;是所有樣本的預測值的平均值.

2 結果與討論

2.1 不同CE 特征模型的預測結果

本研究構建了8 個基于不同組合方式的CE 機器學習特征模型， 其中包括4 個sCE 特征模型(CE-1， CE-3， CE-6， CE-7)和4 個dCE 特征模型(CE-2， CE-4， CE-5， CE-8).這些CE 特征模型對形成能預測的相關系數R2和RMSE 如表2 所示， 其中環境原子選取從中心原子的第1 近鄰到第10 近鄰的原子， 收斂結果以粗體顯示.

表2 8 種CE 特征模型預測的形成能測試集的相關系數R2 和RMSETable 2 Correlation coefficient R2 and RMSE of the predictions of formation energies in the testing data sets by the eight CE feature models

通過對4 個sCE 特征模型的預測性能進行比較發現: ①標量特征模型CE-1[TC+TE]預測的RMSE 在第2 近鄰處收斂于0.252 eV/atom， 而CE-3[OC+OE]預測的RMSE 在第4 近鄰處收斂于0.166 eV/atom， 這表明T 位點標量特征模型的收斂速度比O 位點標量特征模型快， 但是后者的預測精度較高; ②二維矢量特征模型CE-6[TC， TE]預測的RMSE 收斂于第2 近鄰， 而CE-7[OC， OE]預測的RMSE 在第4 近鄰處收斂， 二者預測的RMSE 均約為0.145 eV/atom.可見， T 位點矢量特征模型的收斂速度仍然比O 位點矢量特征模型快， 但二者達到了相似的精度，且均優于標量特征模型.綜合考慮收斂性和準確性， CE-6 是最佳的sCE 特征模型.

通過對4 個dCE 特征模型的預測性能進行比較發現: ①標量特征模型CE-5[TC+TE+OC+OE]預測的RMSE 在第2 近鄰處收斂于0.164 eV/atom; ②矢量特征模型CE-2， CE-4 和CE-8 預測的RMSE 均收斂于約0.145 eV/atom， 優于標量特征模型.二維矢量特征模型CE-2[TC+TE， OC+OE]預測的RMSE 收斂于第1 近鄰， 而CE-4[TC+OC， TE+OE]預測的RMSE 收斂于第2 近鄰.這表明將中心原子與環境原子分離的收斂速度快于將T 位點原子與O 位點原子分離的速度.四維矢量特征模型CE-8[TC， TE， OC， OE]預測的RMSE 收斂于第1 近鄰， 類似于CE-2[TC+TE， OC+OE].因此， 根據收斂性和準確性結果， CE-2[TC+TE，OC+OE]和CE-8[TC， TE， OC， OE]在所有8 個CE 特征模型中具有相似的最佳性能.此外， 四維矢量特征模型CE-8[TC， TE， OC， OE]同時區分了中心原子、環境原子以及T 位點、O 位點，其優點是可以避免因交換T 位點和O 位點異種元素而產生的特征簡并， 但在二維矢量特征模型CE-2[TC+TE， OC+OE]中則會出現這種特征簡并情況.

本研究中的四維和二維矢量CE 特征模型具有相似的預測精度.這可能是因為: 一方面很多T 位點和O 位點占據元素的穩定性， 具有近似交換不變的對稱性; 如果這種位點占據穩定性的交換對稱被破壞， 則四維矢量特征模型應該比二維矢量特征模型具有更準確的描述能力和預測準確性; 另一方面， 如果訓練數據數量有限， 比二維矢量特征模型更復雜的四維矢量特征模型可能更容易導致過擬合.因此， 在CE 特征模型的實際應用中應根據研究結構的復雜性與可用數據集的大小折衷選擇矢量特征模型的維數.如果結構復雜并且數據集很大， 高維CE 特征模型在準確性和可遷移性方面的表現會更好; 反之， 低維CE 特征模型就能達到足夠高的預測精度， 且避免了過擬合.

所有8 種CE 特征模型對晶格常數預測的收斂行為類似于形成能， 結果如表3 所示， 其中環境原子選取從中心原子的第1 近鄰到第10 近鄰的原子， 收斂結果以粗體顯示.但是， 對于晶格常數大于0.7 nm 的結構的預測精度不如其他結構.這可能是由于晶格常數大于0.7 nm 的結構數量較少， 從而導致機器學習模型不能對這類結構進行很好的學習.綜合考慮預測精度和收斂速度， CE-8 特征模型預測的RMSE 在第1 近鄰處達到了0.087°A， 表現最佳.

表3 8 種CE 特征模型預測的晶格常數測試集的相關系數R2 和RMSETable 3 Correlation coefficient R2 and RMSE of the prediction of lattice constants in the testing data sets by the eight CE feature models

2.2 預測結果的統計直方圖

綜合考慮晶格常數預測結果的精度和收斂速度， CE-8 特征模型在最近鄰處的預測結果最好.圖3 和4 是對5 329 種尖晶石氧化物結構形成能和晶格常數的DFT 計算結果的統計分析.5 329 種尖晶石氧化物結構中有141 種是已被實驗合成的結構[19].圖5 是CE-8 特征模型在最近鄰處機器學習預測結果的統計， 其中141 個已知實驗結構的形成能以紅色表示.本研究對這些實驗結構進行了DFT 計算， 發現其形成能在-2.999～-0.633 eV/atom 之間， 因此可將研究的尖晶石氧化物結構分為3 類: ①形成能大于-2.999 eV/atom 的更穩定(more stable，MS)結構; ②形成能介于-2.999～-0.633 eV/atom 之間的穩定(stable， S)結構; ③形成能小于-0.633 eV/atom 的較不穩定(less stable， LS)結構.如圖5(a)所示， 本研究中的尖晶石氧化物結構可劃分成361 個MS 結構、4 347 個S 結構和466 個LS 結構.此外， 還有155 個結構幾何弛豫不收斂， 主要是含有活性金屬Cs 和重稀土元素， 如Eu， Dy， Er 等.

圖3 8 種CE 特征模型預測的形成能與DFT 計算結果的比較Fig.3 Comparisons of the formation energies between the DFT calculations and the predictions of the eight CE feature models

形成能可衡量物質的結合強度， 能作為判斷物質是否熱力學穩定的重要判據之一.本研究通過DFT 計算預測得到了形成能最大、熱力學最穩定的前20 個尖晶石氧化物， 按熱力學穩定性降低順序分別是YbLu2O4，YbTm2O4，YbEr2O4，YbSc2O4，CaLu2O4，YbHo2O4，YbY2O4，CaTm2O4， YbDy2O4， SrLu2O4， CaEr2O4， YbTb2O4， HfYb2O4， CaHo2O4， SrTm2O4，ThYb2O4， CaSc2O4， EuHo2O4， YbGd2O4和CaY2O4.這些通過DFT 計算預測得到的最穩定的尖晶石氧化物結構的形成能范圍是-3.57～-3.38 eV/atom， 遠大于本研究實驗結構的形成能范圍(-2.999～-0.633 eV/atom)， 所以這些預測得到的穩定尖晶石氧化物結構更有可能被實驗合成.本研究的預測結果為將來實驗合成新型的、穩定的尖晶石氧化物結構指明了方向.

本研究預測的尖晶石氧化物的晶格常數(見圖5(b))也呈正態分布， 絕大多數結構的晶格常數處于0.6～0.7 nm 之間.

圖5 CE-8 特征模型預測的形成能和晶格常數的統計分布直方圖Fig.5 Statistical distribution histogram of the formation energies and lattice constants predicted by CE-8 feature model

2.3 機器學習預測的形成能和晶格常數熱圖

機器學習預測的形成能分布如圖6 所示， 其中橫軸和縱軸分別按照T 位點和O 位點組成元素的原子序數排序.根據形成能范圍定義的MS， S 和LS 結構在熱圖中分別大致對應藍色、綠色和橘紅色區域(白色區域是計算不收斂結構).這3 種區域呈現出局部聚集的塊狀， 而非稀疏的點樣分布， 表明具有相近原子序數的結構元素通常形成具有相似穩定性的尖晶石氧化物.此外， 除少數情況外， 大多數聚集的塊狀呈對角對稱分布.這表明將T 位點或O 位點元素替換為鄰近的元素可以形成穩定性相似的結構， 即穩定性的位點互換不變對稱性.這也解釋了為什么會存在很多穩定的反尖晶石和混合尖晶石氧化物.

圖4 8 種CE 特征模型預測的晶格常數與DFT 計算結果的比較Fig.4 Comparisons of the lattice constants between the DFT calculations and the predictions of the eight CE feature models

為便于指導實驗合成新型穩定的尖晶石氧化物， 本研究分別針對T 位點和O 位點將73 個置換元素定性地劃分為與穩定性有關的“好”和“壞”元素(見圖2).如圖6 所示: T 位點和O 位點的穩定性“好”元素組合可形成更穩定的MS 結構， 顯示為圖6 中的藍色區域; T 位點和O 位點的穩定性“壞”元素組合形成了更不穩定的LS 結構， 顯示為圖6 中的橘紅色區域; 如果在T 位點和O 位點分別存在穩定性“好”元素和“壞”元素， 則形成一般穩定的S 結構， 顯示為圖6 中的綠色區域，其穩定性與已知的實驗結構相當.穩定性“好”和“壞”元素都呈現出T/O 位點的對稱穩定性， 即大多數元素對于T 位點和O 位點同時都是穩定性“好”元素或“壞”元素(見圖2 和6).但也有一些少數例外: 例如， Ru 是T 位點的穩定性“好”元素， 但卻是O 位點的穩定性“壞”元素; Sc， Hf， Pa， U 和Np 是O 位點的穩定性“好”元素， 但卻是T 位點的穩定性“壞”元素.

圖6 機器學習預測的形成能熱圖Fig.6 Heat map of the formation energies predicted by ML

根據A， B 陽離子的平均離子半徑， 對預測的尖晶石晶格常數進行約化， 即將晶格常數除以(RA+RB)/2， 其中RA和RB分別表示A， B 陽離子的離子半徑.圖7 為根據A， B 陽離子的平均離子半徑約化后的晶格常數熱圖.可以看到: 約化后的晶格常數熱圖顯示出分塊分布的特性， 即原子序數相近的元素形成的尖晶石氧化物結構具有近似的約化晶格常數; 約化后的晶格常數顯示出沿對角線的對稱性， 表明將T 位點或O 位點元素替換為鄰近的元素可以具有相近的約化晶格常數的結構.相比于同周期的其他元素， 含堿金屬元素Li， Na， K， Rb， Cs 的尖晶石結構通常具有相對較小的約化晶格常數.這可能是由于堿金屬元素相比同周期元素擁有更大的離子半徑.此外， 當A， B 陽離子的原子序數均較小時， 其晶體結構比原子序數均較大時具有明顯更大的約化晶格常數.這可能是由于原子序數均較小時的元素具有較小的電負性， 更容易向氧轉移電子， 從而導致帶更多負電的氧離子間排斥力增加， 約化晶格常數增加.

圖7 根據A， B 陽離子的平均離子半徑約化后的晶格常數熱圖Fig.7 Heat map of the lattice constants after reduction of the average ion radius of A and B cations

2.4 計算和機器學習預測結果與已被實驗合成的尖晶石氧化物的性質比較

為了進一步獨立測試驗證本研究結果的準確性， 將文獻報道的已被實驗合成的141 個尖晶石氧化物結構的晶格常數實驗值分別與DFT 計算值和隨機森林機器學習預測值進行比較， 結果如圖8 所示.可見， 機器學習預測點陣常數的RMSE 為0.04 °A， 比DFT 計算的RMSE(0.09 °A)更低， 進一步證實了機器學習模型的預測有效性.表4 和5 是部分已知尖晶石氧化物結構的形成能和帶隙的實驗值、DFT 計算值以及機器學習預測值.DFT 計算值與部分實驗值的偏差可能來源于交換相關密度泛函的近似.由于將DFT 計算值作為目標值進行機器學習訓練， 所以機器學習和DFT 計算結果比較一致.

表4 部分已知尖晶石氧化物結構的形成能的實驗值、DFT 計算值和機器學習預測值Table 4 Experimental values， DFT calculation values， and ML prediction values of the formation energies of some known spinel oxides (eV·atom-1)

圖8 已被實驗合成的141 個尖晶石氧化物結構的晶格常數實驗值與DFT 計算值和機器學習預測值的比較Fig.8 Comparisons of the experimental lattice constants of 141 known spinel oxide structures with the results calculated by DFT and predicted by ML

表5 部分已知尖晶石氧化物帶隙的實驗值、DFT 計算值和機器學習預測值Table 5 Experimental values， DFT calculation values， and ML prediction values of the band gap of some known spinel oxides eV

3 結束語

本研究通過第一性原理計算和機器學習方法， 系統研究了包含73 種A/B 構成元素的5 329 種AB2O4正尖晶石氧化物.根據141 個已知的尖晶石氧化物結構的形成能范圍(-2.999～-0.633 eV/atom)， 將預測結構分為361 個更穩定(MS)結構、4 347 個穩定(S)結構和466 個較不穩定(LS)結構.為指導實驗成分設計， 提出了與尖晶石氧化物的穩定性定性相關的“好”和“壞”的構成元素.常見的穩定性“好”元素是Ⅱ， Ⅲ， Ⅳ族元素和稀土元素， 大部分具有對稱的占位穩定性.含IA 族堿金屬元素的尖晶石氧化物結構通常具有較大的晶格常數.本研究預測的穩定性尖晶石氧化物結構對進一步實驗合成新型尖晶石氧化物具有指導意義.

本研究構建了包含成分和結構信息的CE 機器學習特征模型.該模型能夠高效預測尖晶石氧化物的能量和結構特性， 并獲得與第一性原理方法相當的精度.通過考察不同的特征構建方式， 發現雙中心CE(dCE)模型的預測性能通常優于單中心CE(sCE)模型， 并且應綜合考慮結構的復雜性和可用數據的數量來選擇復合特征向量的維數.如果結構復雜、對稱性很低并且具有大量數據， dCE 或多中心CE 特征模型的矢量形式將比sCE 模型更加準確可靠.本研究中的CE 特征模型可推廣用于其他凝聚態材料能量和結構的機器學習預測.