氮化硅微波高溫介電函數深度學習分子動力學模擬

李志強 譚曉瑜 段忻磊 張敬義 楊家躍?

1) (山東大學前沿交叉科學青島研究院,光-熱輻射研究中心,青島 266237)

2) (山東大學能源與動力工程學院,濟南 250061)

3) (航天材料及工藝研究所,先進功能復合材料技術重點實驗室,北京 100076)

氮化硅(β-Si3N4)是當下最具應用前景的熱透波材料,其基礎物性高溫介電函數的精準測量對加快氮化硅基熱透波材料的設計,解決高超聲速飛行器“黑障”問題具有重要意義.本文以第一性原理數據為基本輸入,基于深度神經網絡訓練得到深度學習勢,然后運用深度學習分子動力學方法計算了氮化硅高溫微波介電函數.與傳統經驗勢相比,深度學習勢的計算結果與實驗結果在數量級上保持一致;同時發現,深度學習分子動力學在計算速度方面表現優異.此外,建立了氮化硅弛豫時間溫度依變性的物理模型,揭示了弛豫時間溫度依變性規律.本研究通過實現大規模高精度的分子動力學模擬計算了氮化硅高溫微波介電函數,為推動氮化硅基材料在高溫熱透波領域的應用提供了基礎數據支撐.

1 引言

飛行器高超聲速飛行時經受劇烈的氣動加熱,表面溫度高達700—3000 K,天線罩會出現“黑障”(blackout)[1?3]現象.為避免“黑障”現象影響通信質量,要求熱透波材料在高溫狀態下有極低的介電參數和介電損耗,且在很寬的溫度區間內變化幅度小.氮化硅具有優異的高溫力學、抗熱震性能以及良好的介電特性,已成為一種具有重要應用前景的高溫透波材料.然而,由于實驗條件的限制,很難在極端高溫條件下測量氮化硅的介電性能.另外,飛行器價值不菲,如僅憑實驗測量來收集數據參數,成本較高,并且由于高速飛行器回收困難,其原始實驗數據的準確性很難保證.鑒于以上技術難點,迫切需要探索一種精確計算氮化硅高溫介電參數的理論方法.

當前氮化硅介電特性的研究主要集中在采用實驗方法結合電介質理論建立理論模型對高溫下氮化硅介電特性進行推算[4,5].理論模型中的關鍵參數依賴實驗測量或經驗假設,難以適用于計算熱透波材料的高溫微波介電函數.分子動力學(molecular dynamics,MD)模擬可從原子尺度準確描述材料的結構演變過程和計算相關物性,是獲得電介質理論關鍵參數的重要手段.研究表明[6?8],基于線性響應理論可將分子動力學模擬計算得到的偶極矩數據用于計算材料介電常數.Afify 等[9]和Cardona 等[10]利用分子動力學模擬研究了水、乙醇、乙二醇以及乙醇胺的微波介電特性,發現極化經驗力場更適合用于描述物質的介電特性.由此可知力場是保證MD 模擬準確的重要前提,經驗力場是由實驗數據和第一性原理數據擬合得到,顯然不適用于描述氮化硅高溫微波條件下原子間相互作用.機器學習勢[11,12]結合AIMD (ab initiomolecular dynamics)方法精度高和MD 方法計算效率快的優點,從求解電子薛定諤方程出發構造原子間相互作用勢而無需引入經驗模型,為計算氮化硅高溫微波介電函數指明了研究方向.

目前常見的機器學習勢有高斯近似勢[13]、高維神經網絡勢[14]、矩張量勢[15]和深度學習勢 (deep learning potential,DLP)[16]等方法.得益于圖形處理單元 (graphics processing unit,GPU)的加速計算,DLP 在計算效率和精度方面有更出色的表現.Chen 和Li[17]訓練了DLP 用于研究300 K 溫度下碳化硅的介電性質,計算結果與實驗測量數據吻合良好.然而,應用DLP 研究氮化硅高溫介電性質的工作卻鮮有報道.

本文利用第一性原理數據訓練了DLP,并應用深度學習分子動力學 (deep potential molecular dynamics,DPMD)方法模擬了氮化硅的極化弛豫過程.DLP 由不同條件下氮化硅結構訓練,其精度通過比較第一性原理計算結果與預測結果來確定.基于訓練的DLP,本文計算了氮化硅的微波高溫介電函數,計算結果與實驗結果高度吻合.此外,建立了氮化硅弛豫時間溫度依變性的數學模型,理解和掌握了弛豫時間隨溫度的變化規律.

2 模型和方法

本文采用融合了電子結構、機器學習和分子動力學3 種手段的多尺度方法研究了氮化硅的微波高溫介電性質,技術路線如圖1 所示.首先,進行小規模的AIMD 計算用于采集指定熱力學條件下的原子結構信息;其次,以選取原子結構的能量、受力、體系尺寸以及原子位置信息為輸入,經深度神經網絡 (deep neural network,DNN) 進行訓練擬合后得到深度學習勢;最后,基于深度學習勢進行大規模高精度的分子動力學模擬.

圖1 深度學習勢訓練流程: 主要包含AIMD 計算采樣、深度神經網絡訓練以及深度學習分子動力學模擬三部分Fig.1.DLP training process: AIMD computational sampling,DNN training,and DPMD simulation.

2.1 訓練集獲取

深度學習勢函數的精度取決于訓練集的合理性和豐富性.本文使用CP2K[18]軟件來生成初始訓練構型,初始訓練結構的體系原子數目為140,如圖2(a)所示.具體過程分為以下3 步: 首先,利用Broyden-Fletcher-Goldfarb-Shanno (BFGS) 算法對β-Si3N4初始結構進行結構優化;其次,整個體系在正則系綜(canonical ensemble,NVT)下通過速度標定法進行熱力學弛豫平衡,溫度設置為300—1000 K (溫度間隔為100 K),時間步長為0.5 fs;最后,每個溫度條件下各選取一段0.5 ps 的平衡軌跡并記錄每個時間步下的體系能量、體系尺寸、原子坐標以及受力信息.混合高斯平面波(Gaussian and plane wave,GPW)[19,20]的密度泛函理論用于評估能量與力: 計算基組選用DZVP-GTHPADE,贗勢則分別選用GTH-PADE-q4 (Si)和GTH-PADE-q5 (N),截斷能設置為400 Ry,k點網格設置為1×1×1.樣本空間的豐富性是決定勢函數訓練結果好壞的先決條件,因此對AIMD平衡軌跡中的構型進行篩選并剔除重復構型,最終共計選出871 個結構用于構建深度學習勢函數,隨機選取其中88 個(占比10%)結構作為驗證集,剩余783 個(占比90%)結構用于勢函數的訓練.

圖2 β-Si3N4 的晶體結構 (a) 用于第一性原理分子動力學計算的體系(原子個數為140 個);(b) 用于深度學習分子動力學模擬的體系(原子個數為2200 個),其中綠色原子代表硅原子,紅色原子代表氮原子Fig.2.Crystal structure of β-Si3N4: (a) The system for first-principles molecular dynamics calculations (140 atoms);(b) the system for deep potential molecular dynamics simulations (2200 atoms).The green atoms represent silicon atoms and the red atoms for nitrogen atoms.

2.2 深度學習勢函數訓練

本文采用DeepMD-kit[21]軟件包訓練深度學習勢,該軟件基于當下流行的深度學習框架Tensorflow[22]進行二次開發,能夠實現多體勢訓練并保持計算精度幾乎不損失.依據原子神經網絡勢的理論框架[12,23]可知,對于包含N個原子的體系,體系總能量可歸結為每個原子的能量之和:

其中每個原子的能量Ei由該原子及其截斷半徑范圍內的近鄰原子位置所決定:

其中,Ri表示原子i的原子位置,Rj代表原子i的近鄰原子j的原子位置,代表原子i的截斷半徑范圍內近鄰原子列表的索引集.

其中ps,pf和pξ分別代表能量項、受力項及維里項的可調前置因子;Δ 代表訓練集數據與深度學習勢預測值間的偏差;N代表原子數目;ε代表單個原子的能量;Fi代表第i個原子的受力信息;ξ代表維里張量.由于訓練過程中沒有包括維里值數據,將pξ設為0.對于損失函數中的能量和力誤差,可調預因子的起始值分別為0.02 和1000,這意味著訓練過程在初始階段主要由力擬合控制.學習率從最初的1.0×10–3以指數形式下降到3.5×10–8,衰減率和衰減步長分別為0.95 和5000,訓練步數設置為400 萬步.

2.3 深度學習分子動力學模擬

借助DeepMD-kit 與LAMMPS[26](large-scale atomic/molecular massively parallel simulator)之間的軟件接口,可以直接實現基于深度學習勢的分子動力學模擬.本文使用深度學習勢對含2200 個原子的β-Si3N4體系進行分子動力學模擬計算,具體實施過程如下: 體系中原子依據高斯分布進行速度初始化,時間步長設置為0.5 fs.首先,整個β-Si3N4體系在對應溫度下進行熱力學弛豫,弛豫時間為1 ns;其次,待體系達到熱力學穩定狀態后,開始每2000 步收集一次偶極矩數據用于后續計算,采集時長為1 ns.上述弛豫過程和數據采集過程均在NVT 系綜完成,溫度控制由Nosé-Hoover[27,28]熱浴實現.為了比較經驗勢函數和深度學習勢的差異,同樣使用Tersoff[29,30]經驗勢進行了上述模擬.

3 結果與討論

3.1 勢函數訓練及結果驗證

為了驗證訓練結果的準確性,本文以均方根誤差(root mean square error,RMSE)來衡量深度學習勢的預測結果與AIMD 計算結果之間的差異.圖3(a)為DLP 預測能量值與AIMD 計算值的對比,其中數據點沿對角線分布表明深度學習勢的預測能力較好,能量的RMSE 為0.00550 meV/atom;圖3(b)—(d)為x,y,z三個方向上DLP 原子受力預測值與AIMD 計算值對比,圖中紅色實線代表預測結果與計算結果相等,數據點離直線越近說明計算誤差越小,RMSE 結果為7.800 meV/?.由圖3 可知,體系能量和原子力的數據點絕大多數都落在直線上,說明深度學習勢的預測能力較好,預測結果已經達到了AIMD 計算精度.

圖3 (a) 第一性原理計算體系能量和深度學習勢計算體系能量對比關系圖;(b)—(d) 第一性原理計算原子受力和深度學習勢計算原子受力對比關系圖.其中圖中直線代表y=xFig.3.(a) Comparison between the energy calculated by first-principles and that by the deep learning potential;(b)–(d) comparison between the forces on the atoms calculated by first-principles and those by the deep learning potential.The straight line in the figure represents y=x.

為進一步評估深度學習勢相對于AIMD 計算的準確性,本文分別使用AIMD 方法和DPMD 方法計算了體系的徑向分布函數進行對比,來驗證深度學習勢在描述結構演化過程中的精度和準確性.在相同溫度條件下以AIMD 初始結構進行深度學習分子動力學模擬,計算體系原子數目均為140.兩種模擬得到的Si-Si,Si-N,N-N 的徑向分布函數如圖4 所示,可以明顯看到深度學習勢準確地復現了AIMD 徑向分布函數峰的形狀和位置,表明深度學習勢函數可以很好地描述體系的動力學性質.

圖4 (a)—(c) 300 K 溫度下DPMD 與AIMD 計算徑向分布函數對比圖;(d)—(f) 1000 K 溫度下DPMD 與AIMD 模擬徑向分布函數對比圖Fig.4.(a)?(c) Comparison of radial distribution function between DPMD and AIMD simulations at 300 K;(d)?(f) comparison of radial distribution function between DPMD and AIMD simulations at 1000 K.

3.2 弛豫時間及微波介電函數計算

本文依據德拜方程的修正式柯爾-柯爾[31](Cole-Cole)公式來計算β-Si3N4的微波段介電函數:

其中,ε∞為光頻介電常數,εs為靜態介電常數,ω為電磁波的角頻率,τ為弛豫時間,α為弛豫修正因子.對應的介電函數實部和虛部分別表示為

式中,ε′為介電常數的實部,ε′′為介電常數的虛部.據(5)—(7)式可知,計算氮化硅材料頻率相關介電函數需要確定弛豫時間、光頻介電常數、靜態介電常數以及弛豫修正因子等參數.由于高頻與靜態介電常數受溫度影響較小,故采用0 K 溫度下氮化硅的光頻和靜態介電常數用于計算,其值分別為4.27 和8.16[32].弛豫修正因子通常用于修正理論值與實驗值間的差異,本文主要采用模擬計算的方法計算介電函數,因此將弛豫修正因子設置為0.弛豫時間能通過求解偶極矩自相關函數得到[17],偶極矩自相關函數計算如下:

式中M(t)代表某時刻體系的總偶極矩,M(0)為初始時刻體系的總偶極矩.依據德拜模型,自相關函數通常以指數遞減形式呈現,偶極矩自相關函數與弛豫時間存在如下關系:

其中t為時間.

圖5 為偶極矩自相關函數隨時間變化曲線,將偶極矩自相關函數擬合成(9)式中的指數衰減形式,以確定弛豫時間τ.擬合過程如下: 首先通過分析得出偶極矩自相關函數存在局部最大值,即圖5(a)中曲線的所有峰值.然后將得到的極值點以對數形式進行擬合,擬合結果如圖5(b)所示,所得直線的斜率的倒數即為弛豫時間.偶極矩數據是按照一定頻率進行采樣的,部分數據并不能完全落在自相關函數的極值點處.因此,本文選擇了數值大于其最近鄰值的點作為極值的近似,這也會導致弛豫時間計算時中擬合數據與原始數據不能完全重合.通過上述方法計算了300—1000 K 溫度范圍內的弛豫時間,具體結果匯總在表1 中.結果發現,利用深度學習勢與經驗勢得到的弛豫時間在數量級上相當,且隨著溫度的升高,弛豫時間均呈現下降的趨勢,這是由于弛豫時間遵循阿倫尼烏斯(Arrhenius)定律,隨溫度升高呈下降的趨勢[33].值得注意的是,Tersoff 勢函數計算得到的弛豫時間絕大多數情況下大于DLP 的計算結果,但是結果的波動性較大.原因可能是Tersoff 勢最多描述三體相互作用,其勢函數形式只是粗略的擬合近似.深度神經網絡在理論上可以用來逼近任何高維函數,具有表征多體勢的優點,因此由DLP 計算得到的結果更加準確,所表現的物理圖像更具說服力.

圖5 運用深度學習分子動力學計算的700 K 溫度下氮化硅偶極矩自相關函數 (a) 偶極矩自相關函數隨時間的變化;(b) 偶極矩自相關極值點對數值隨時間的變化Fig.5.Silicon nitride dipole moment autocorrelation function at 700 K calculated by the deep learning potential: (a) Dipole moment autocorrelation function versus time;(b) dipole moment autocorrelation polar point logarithm versus time.

根據(5)—(7)式和表1 中弛豫時間即可求得氮化硅的微波段 (0—100 GHz)介電函數.圖6 展示了β-Si3N4在300—1000 K 溫度范圍內頻率相關介電函數的實部和虛部.其中文獻[34]結果為300—800 K 溫度時在頻率為8.2 GHz,10 GHz,12.4 GHz 時的介電函數實部.300 K 溫度下,DLP和Tersoff 勢在頻率為8.2 GHz 時結果與實驗值相差不大,但隨著頻率的升高,DLP 的結果與實驗值更加接近.文獻[35]為常溫下氮化硅摻雜多壁碳納米管后介電函數實部變化,隨著頻率的升高介電函數實部呈現下降的趨勢,這與本文計算結果的變化趨勢是相符的.因此,隨著頻率增大,介電函數實部呈現逐漸降低的趨勢,DLP 的計算結果均大于Tersoff 勢結果,但與實驗值在數量級上更加吻合.從圖6(b)給出的微波介電函數虛部來看,300 K溫度下DLP 計算的介電函數虛部相比Tersoff 勢更接近于實驗值,表明了DLP 可以更好地描述β-Si3N4的高溫微波介電特性.

表1 不同溫度條件下的弛豫時間Table 1.Relaxation time under different temperatures by DLP and Tersoff potential.

圖6 (a) β-Si3N4 在不同溫度時的頻率相關介電函數實部;(b) β-Si3N4 在不同溫度時的頻率相關介電函數虛部.其中文獻[34]僅給出了實部值,未給出虛部值.文獻[35]給出了8—12 GHz 頻率范圍內的介電函數實驗測量值Fig.6.(a) Real part and (b) imaginary part of frequency-dependent dielectric function of β-Si3N4 at varying temperatures.Note that only the real part is given in Ref.[34].The values of dielectric function in the frequency range of 8–12 GHz are given in Ref.[35].

對于Tersoff 勢來說,其構造過程是先選擇合適的數學表達式,然后通過實驗數據擬合確定函數中的參數,因此經驗勢只能描述已知的性質.因此,利用Tersoff 勢來計算高溫氮化硅的介電特性還有待商榷.相比之下,基于第一性原理計算求解電子結構所得到結果是值得信賴的,而深度學習勢是基于第一性原理數據進行訓練擬合,因此深度學習勢計算得到的氮化硅介電函數結果更加合理.

3.3 弛豫時間溫度依變性模型

由德拜模型可知,弛豫時間是計算微波段介電函數的關鍵參數,掌握弛豫時間隨溫度變化規律是準確計算高溫微波介電函數的重要前提.本文基于深度學習勢弛豫時間計算結果,擬合得到了弛豫時間隨溫度變化曲線:

由圖7 可知弛豫時間的對數隨溫度變化呈現線性下降的趨勢且直線相關性較好,因此可判定此模型在對應溫度范圍內的有效性.結合(10)式和圖7 可知,隨著溫度的升高,弛豫時間呈現指數下降的趨勢.當溫度高于1700 K,弛豫時間趨近于零.這是因為高溫情況下極化損耗占比較少,需要考慮電子電導以及等離子的損耗貢獻,這與文獻[36]的報道結果一致.

圖7 弛豫時間隨溫度變化曲線,其中藍色直線代表弛豫時間溫度依變性模型,紅色散點代表不同溫度下弛豫時間變化曲線Fig.7.Relaxation time variation curves with temperature.The blue straight line corresponds to relaxation time temperature dependence model,and the red scattered points represent relaxation time variation curves at different temperatures.

3.4 計算效率測試

DPMD 方法實現了具有第一性原理高精度的大規模分子動力學模擬,在計算效率上也有了很大的提升.由于AIMD 方法很難實現數千原子體系的計算求解,故而此處選用原子數為140 的氮化硅體系分別測試了DPMD 模擬與AIMD 模擬的計算效率.本文采用了48 個Intel Xeon Platinum 9242 CPU 核心通過并行計算方式實現AIMD 模擬計算.DeepMD-kit 軟件有CPU 和GPU 兩種版本,本文選用后者并在Nvidia RTX 3080 進行了模擬計算,計算結果如圖8 所示.可以看出,得益于GPU 加速,DPMD 模擬的計算速度大幅度提升,表明了DLP 在GPU 加速的條件下既能保持AIMD模擬的計算精度,又能在計算速度上有優異的表現.

圖8 DPMD (Nvidia RTX 3080 GPU 計算)與AIMD(48 個Intel Xeon Platinum 9242 CPU 計算)計算速度結果Fig.8.Computational speed of DPMD (running with a Nvidia RTX 3080 GPU) and AIMD calculations (running with 48 Intel Xeon Platinum 9242 CPU cores).

4 結論

本文基于DNN 對高溫氮化硅分子晶體結構數據集進行訓練得到深度學習勢,并應用DPMD方法結合電介質理論研究了β-Si3N4的高溫微波介電函數.相較于第一性原理計算結果,深度學習勢所計算原子受力的RMSE 為7.800 meV/?,體系能量的RMSE 為 0.00550 meV/atom,均能很好地預測正確的結果.由于深度學習勢良好的訓練結果,DPMD 在模擬氮化硅體系結構演化過程時的徑向分布函數與AIMD 結果完全一致.同時發現,深度學習勢計算得到的氮化硅高溫微波介電函數與實驗結果在數量級上保持一致.此外,基于深度學習勢弛豫時間的計算結果,建立了弛豫時間隨溫度變化的數學模型,探明了弛豫時間的溫度依變性規律.本文研究結果成功預測了氮化硅高溫微波介電函數,為推動其在高溫熱透波領域的應用提供了基礎數據支撐和科學依據.