基于物理模型的BaZrO3鈣鈦礦機器學習力場

趙 亮, 牛宏偉, 荊宇航

（1.中國航發北京航空材料研究院，北京 100095；2.哈爾濱工業大學 航天學院，哈爾濱 150001）

高性能航空發動機往往要承受更高的渦輪進口溫度。目前，推重比10一級的航空發動機渦輪前溫度已經超過1800 K，推重比的提高將進一步提升渦輪前溫度。鎳基高溫合金作為當前渦輪葉片的主要材料，其耐高溫能力只有1373 K，因此為了保證長期可靠工作，渦輪葉片普遍由鎳基單晶基體、熱障涂層（TBC）以及復雜氣冷結構組成。其中，熱障涂層在提高渦輪前溫度和延長葉片使用壽命方面成效顯著，已成為高性能發動機研制的關鍵技術之一。20世紀，美、英、法、日、俄等國家都在積極研究熱障涂層的設計與制備技術，并大量應用于航空發動機熱端部件。進入21世紀，熱障涂層技術逐漸成熟并得到廣泛應用。俄羅斯的蘇-30和蘇-35戰斗機的噴管都使用了熱障涂層材料，而美國目前幾乎所有的軍用和民用航空發動機都采用了熱障涂層技術[1-4]。

鈣鈦礦[5-8]結構化合物具有高熔點、低熱導率以及高膨脹系數等特點，是熱障涂層陶瓷的備選材料之一，典型的包括BaZrO3（BZO）等。很多學者對BZO進行了大量的密度泛函理論（DFT）計算研究，例如其結構、熱力學和質子傳輸方面，特別是在摻雜狀態下[9]，因為摻雜物會引起結構畸變并影響質子傳輸能力。BZO基態結構的中子散射實驗和第一原理計算[10]揭示了不穩定的聲子模式是一種過渡現象，這與ZrO6八面體的旋轉有關，這些依賴于交換相關函數的計算表明了BZO的結構是穩定立方鈣鈦礦氧化物。BZO作為一種脆性材料，它的力學行為受到其自身固有微觀結構的強烈影響[11-12]，如晶界、預先存在的微裂紋、晶粒和孔隙。這些微結構的取向、尺寸、形狀和空間分布控制著材料的開裂、延伸、相互作用和最終失效[13-14]。因此，這些微觀結構的變化可以改變裂紋的發生和擴展行為，導致結構最終的失效，從而改變材料中韌性和斷裂強度的分散性。然而，為了更深入地理解BZO力學行為的本質，需要原子尺度模擬來獲取必要的微觀信息[15]。DFT[16-17]計算通常用于研究小系統的靜態特性，而分子動力學（MD）模擬[18-19]已被廣泛用于研究更大系統的動力學特性。力場對分子動力學計算結果起決定性作用，力場的精確度決定了分子動力學模擬的精確度。BZO體系目前已考慮兩體短程作用和庫侖靜電作用的經驗力場，但力場的計算精度有待提高。

廣義來講，MD力場分為基于物理模型的“純物理”的經驗力場[20-24]和基于機器學習算法的“純數學”的機器學習力場[25-30]。在機器學習力場出現之前，經驗力場是分子動力學模擬的主導。針對各種不同的體系和不同的科學問題，研究者已經開發出各種各樣的經驗力場，這些經驗力場凝結了前人對于模擬體系和科學問題的物理理解。經驗力場的優勢是可移植性較好，因為其中包含了對體系的物理描述。但這些經驗力場中的物理假設缺陷在于沒有完全描述體系中的所有物理本質，例如有的經驗力場僅考慮了體系的兩體作用和三體作用，而忽略了體系的其他多體作用，勢能面計算精度較低?！凹償祵W”的機器學習力場不含有任何對體系的物理假設[31-34]，通過大量的擬合參數和靈活的函數形式實現了極強的表示能力，因此勢能面計算精度較高。但機器學習力場最大的缺陷在于可移植性較差：對于和訓練集描述的構象空間相差較大的構象，機器學習力場可能會給出誤差較大的結果。所以，一個自然的想法是，能否將“純物理”的經驗力場和“純數學”的機器學習力場的優點相結合，即實現高精度的勢能面，又實現較好的可移植性。

本研究以BZO體系為例提出了一般的基于物理模型的機器學習力場開發方法。對于經驗力場沒有考慮到的物理作用，采用機器學習的方法進行學習，總的力場就包含了經驗力場和機器學習修正項兩部分。經驗力場考慮了體系主要的物理作用，并且使模型具有了較好的可移植性。機器學習修正項考慮了經驗力場沒有考慮到的其他多體作用，提高了模型的計算精度。這樣就集合“純物理”的經驗力場和“純數學”的機器學習力場的優點，實現高精度高可移植性的力場模型開發，最后從靜態性質、相穩定相和力學性質三方面對基于物理模型的機器學習力場的準確性進行了驗證。

1 機器學習數據集及方法

以BaZrO3體系的靜態性質和力學性質為研究目標，構造了和靜態性質和力學性質密切相關的數據集，包括應變下的BaZrO3構象和沿不同方向單軸拉伸的構象。所有構象均來自于從頭算分子動力學（AIMD）[35]模擬或者基于經驗力場的MD模擬。詳細的數據集如表1所示，該數據集共包含60000個構象。其中，（δ，δ，δ，0，0，0）應變的BaZrO3構象包含了無應變、±1%應變、±3%應變和±5%應變等不同類型，這些構象將用于學習如晶格常數、結合能以及彈性常數等靜態性質；沿不同晶向單軸拉伸的構象將用于學習不同晶向的楊氏模量等力學性質。

本工作使用的經驗力場來源于參考文獻[36]，考慮了兩體短程作用和庫侖靜電作用。經驗力場表達式為：

式中：Uij是原子i和原子j之間的勢能；rij是原子i和原子j之間的距離；qi和qj是原子i和原子j的電荷數；?0為真空介電常數；Aij，ρij和Cij為原子i和原子j之間的Buckingham力場參數。

BaZrO3體系中不同類型原子之間的Buckingham力場參數如表2所示。長程的庫侖靜電作用采用Particle-Particle Particle-Mesh（PPPM）算法[37]求解。

表2 BaZrO3體系中不同類型原子之間的Buckingham力場參數Table 2 Interatomic potentials parameters for BaZrO3

對于數據集中的構象，分別使用DFT和經驗力場[36]計算每個構象的體系勢能和原子受力，兩者的差值被用做機器學習的數據集。最終的力場模型由經驗力場和機器學習力場修正項構成。使用VASP軟件包[38-40]進行DFT計算，計算采用增強投影波方法和Perdew-Burke-Ernzerhof泛函[41]。平面波基組的截斷能設置為400 eV。所有AIMD模擬均使用了2×2×2的Monkhorst-Pack K點網格。圖1給出了構建基于物理模型的機器學習力場的流程。

圖1 構建基于物理模型的機器學習力場的流程Fig. 1 Workflow for constructing a machine learning force field based on a physical model

基于上述數據集，使用DeePMD[42]人工神經網絡框架訓練了機器學習力場。訓練的截斷半徑為0.6 nm，批次大小設置為4，訓練迭代了10萬個批次。起始和終止的學習率分別設置為5×10-3和5×10-8，學習率在訓練過程中指數衰減。最后，使用一個隱藏層結構為240×240×240的深度神經網絡，將局部原子環境映射到原子對體系勢能的貢獻。訓練完成后將獲得基于物理模型的機器學習勢。

2 結果與討論

訓練過程中訓練集體系勢能和原子受力的均方根誤差（RMSE）隨訓練批次的變化如圖2所示，可以看出在訓練結束時體系勢能和原子受力的RMSE均達到收斂狀態。

圖2 訓練過程中訓練集體系勢能和原子受力的RMSE隨訓練批次的變化Fig. 2 Variation of RMSE of potential energy and atomic forces on training set against training batches during training

圖3為數據庫和機器學習力場的體系勢能和原子受力之間的比較。訓練集和測試集的體系勢能RMSE分別為2.09 meV/atom和1.99 meV/atom。對于原子受力，訓練集和測試集的RMSE分別為1182.99 meV/nm和1180.99 meV/nm。機器學習力場的RMSE值與機器學習模型的典型RMSE值一致[43]。

圖3 機器學習力場計算結果和密度泛函理論參考數據的比較 （a）訓練集每個原子平均能量的比較；（b）測試集每個原子平均能量的比較；（c）訓練集上各個原子受力的比較；（d）測試集上各個原子受力的比較；（e）訓練集每個原子的平均能量和各個原子受力的分布；（f）測試集每個原子的平均能量和各個原子受力的分布Fig. 3 Comparisons of machine learning potential results with respect to the references of DFT data （a） comparison of average energy per atom for training set；（b） comparison of average energy per atom for test set；（c）comparison of forces on individual atoms on a training set；（d）comparison of forces on individual atoms on test set；（e） error distributions for energies and forces for training set；（f）error distributions for energies and forces for test set

為了驗證基于物理模型的機器學習力場的準確性，計算了靜態性質，并和DFT理論計算結果進行了比較。這里的靜態性質包括晶格常數a，彈性常數C11、C12和C44。詳細的結果如表3所示。同時，還給出了經驗力場和直接學習訓練集數據得到的純機器學習力場的結果。

表3 基于物理模型的機器學習力場的靜態性質和DFT理論計算結果的比較Table 3 Comparison of static properties based on empirical force field with machine learning corrections and DFT references

從表3可知，本研究的DFT計算結果與其他DFT工作結果一致。另外，無論是基于物理模型的機器學習力場還是純機器學習力場，靜態性質結果都接近DFT參考數據，遠好于經驗力場。最后，對于靜態性質，基于物理模型的機器學習力場獲得了和純機器學習力場接近的結果，說明機器學習修正項可以提高經驗力場計算精度，使其具有類似于純機器學習力場的精度。

合理的力場可以使結構從不穩定轉變為穩定狀態。例如，將不穩定的無序結構轉變為穩定的有序結構，抑或對模擬盒子施加形變后結構可以恢復原狀。為了驗證基于物理模型的機器學習力場的準確性，還進行了如下測試：（1）對BaZrO3體系所有原子施加一個隨機的位移，構造出無序結構，然后進行能量最小化，觀察結構能否轉變為穩定的有序結構恢復原狀；（2）對模擬盒子施加不同的應變，然后進行能量最小化，觀察結構能否恢復原狀。同時，還給出了經驗力場和直接學習訓練集數據得到的純機器學習力場的結果。

分子動力學模擬采用LAMMPS軟件包[14]，三邊均采用周期性邊界條件，能量最小化采用共軛梯度算法。使用徑向分布函數（radial distribution function，RDF）表征結構的有序性。RDF通常指的是給定某個粒子的坐標，其他粒子在空間的分布幾率，因此RDF既可以用來研究物質的有序性，也可以用來描述電子的相關性。

RDF的定義為：

式中：g(r)是徑向分布函數；n(r)是半徑r處，寬度為Δr內的原子總數；ρ是平均原子密度。

對原子數目為625的BaZrO3體系所有原子施加一個大小為0.053 nm（25%的O—Ba鍵長）的隨機位移構造無序結構。圖4為使用基于物理模型的機器學習力場能量最小化過程中結構和RDF的變化（圖中紅色為O原子，藍色為Ba原子，黃色為Zr原子）?？梢园l現隨著能量最小化的進行，結構逐漸從無序轉變為有序，RDF的峰由低變高、由寬變窄。能量最小化結束后，體系恢復為完美的結構。

對原子數目為625的BaZrO3體系所有原子施加了從5%到30%的O—Ba鍵長的隨機位移，分別使用經驗力場、機器學習力場和基于物理模型的機器學習力場進行能量最小化，結果如表4所示?？梢园l現，對于20%及以內的O—Ba鍵長的隨機位移，所有三個模型均可以使無序結構恢復到完美有序結構。但施加25%的O—Ba鍵長的隨機位移時，機器學習力場無法使無序結構恢復到完美有序結構。最后施加30%的O—Ba鍵長的隨機位移時，所有三個模型均無法使無序結構恢復到完美有序結構。由于訓練集多數構象均位于平衡位置附近，所以機器學習力場在描述遠離平衡位置的構象時可能產生較大誤差，比如計算施加很大隨機位移的無序構象時。但經驗力場由于具有物理模型，對于這種遠離平衡位置的構象仍然具有一定的有效性，因此在處理無序構象恢復為完美有序結構時具有優勢?；谖锢砟Ｐ偷臋C器學習力場繼承了經驗力場的這一優勢，在這一場景中展示出了相似的結果，好于純機器學習力場。

表4 基于不同模型的無序結構能量最小化結果Table 4 Energy minimization results for disordered structures based on different models

另外，對原子數目為625的BaZrO3體系模擬盒子施加應變，分別使用經驗力場、純機器學習力場和基于物理模型的機器學習力場進行能量最小化，結果如表5所示?？梢园l現，當同時壓縮模型，并且改變夾角時，純機器學習力場無法使施加應變結構恢復到完美有序結構。此時經驗力場和基于物理模型的機器學習力場仍然有效。這一測試再次說明基于物理模型的機器學習力場繼承了經驗力場中含有物理模型的優勢，能使遠離平衡位置的結構恢復到平衡位置，效果好于純機器學習力場。

表5 基于不同模型的施加應變結構能量最小化結果（其中 為施加應變后的晶格參數，為施加應變前的晶格參數）a，b，c，α，βandγ a0，b0，c0，α0，β0andγ0 a、b、c、α、β和γa0、b0、c0、α0、β0和γ0Table 5 Energy minimization results for structures under strain based on different models（ are lattice parameters after applying strain， are lattice parameters before applying strain）

針對BaZrO3體系，研究了沿不同晶向的線彈性階段的力學性質。表6給出了模擬體系的參數。

表6 BaZrO3模擬體系參數Table 6 Parameters of different BaZrO3 systems

分子動力學模擬采用LAMMPS軟件包[14]。三邊均采用周期性邊界條件模擬宏觀材料在單軸拉伸條件下的線彈性階段的力學性質。拉伸過程采用NPT系綜，控溫算法使用Nose-Hoover熱浴法[45-46]，控壓算法使用Parrinello-Rahman算法[47]。時間步長取0.5 fs。體系首先使用NVT系綜弛豫在300 K下弛豫50 ps，然后每10000步加載0.5%應變。

圖5展示了BaZrO3模型沿不同晶向的線彈性階段應力-應變曲線。表7所示為使用線彈性階段的應力-應變曲線擬合得到的楊氏模量。圖5（a）為模型A在300 K下沿著 [100] 方向的線彈性階段的應力-應變曲線。圖5（b）、（c）和（d）分別為模型B在300 K下沿著方向的線彈性階段的應力-應變曲線?？梢园l現對于所有4個單軸拉伸晶向，機器學習力場和基于物理模型的機器學習力場結果非常接近，而兩者和基于經驗力場的結果存在差異。

圖5 BaZrO3模型沿不同晶向的線彈性階段應力-應變曲線Fig. 5 Stress-strain curves of BaZrO3 in linear elastic region along different crystal directions

表7 BaZrO3模型沿不同晶向的楊氏模量Table 7 Young’s modulus of BaZrO3 along different crystal directions

由表7可知，純機器學習力場和基于物理模型的機器學習力場計算得到的楊氏模量和實驗結果更接近，優于經驗力場計算結果。

3 結論

（1）開發了基于物理模型的機器學習力場。首先構建了數據庫，其中包含了與靜態性質和力學性質相關的構象。然后對于數據集中的構象，分別使用DFT和經驗力場計算了體系勢能和原子受力信息。最后，使用機器學習方法來學習DFT和經驗力場之間的差值。最終的模型由經驗力場和一個機器學習修正項構成。

（2）分別使用經驗力場、純機器學習力場和基于物理模型的機器學習力場計算了BaZrO3的靜態性質，發現純機器學習力場和基于物理模型的機器學習力場的計算結果和DFT一致，遠好于經驗力場。從彈性常數C11、C12和C44的計算來看，純機器學習力場與DFT的計算結果誤差為0.34%、8.75%和10.71%，基于物理模型的機器學習力場與DFT的計算結果誤差為0.34%、2.5%和7.14%。后者的誤差略低于前者，說明機器學習修正項可以提高經驗力場計算精度，使其具有類似于純機器學習力場的精度。

（3）分別使用經驗力場、純機器學習力場和基于物理模型的機器學習力場研究了BaZrO3的相穩定性，發現經驗力場和基于物理模型的機器學習力場在處理無序構象恢復為完美有序結構時具有優勢，基于物理模型的機器學習力場繼承了經驗力場的這一優勢，好于純機器學習力場。

（4）分別使用經驗力場、純機器學習力場和基于物理模型的機器學習力場計算了BaZrO3的四個不同晶向的楊氏模量，發現純機器學習力場和基于物理模型的機器學習力場結果非常接近。純機器學習力場、基于物理模型的機器學習力場和經驗力場計算結果與試驗值的誤差分別為9.22%、1.6%和29.06%，前兩者的誤差要遠低于經驗力場。