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

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

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

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

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

1 引言

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

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

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

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

2 模型和方法

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

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

2.1 訓練集獲取

深度學習勢函數(shù)的精度取決于訓練集的合理性和豐富性.本文使用CP2K[18]軟件來生成初始訓練構(gòu)型,初始訓練結(jié)構(gòu)的體系原子數(shù)目為140,如圖2(a)所示.具體過程分為以下3 步: 首先,利用Broyden-Fletcher-Goldfarb-Shanno (BFGS) 算法對β-Si3N4初始結(jié)構(gòu)進行結(jié)構(gòu)優(yōu)化;其次,整個體系在正則系綜(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點網(wǎng)格設置為1×1×1.樣本空間的豐富性是決定勢函數(shù)訓練結(jié)果好壞的先決條件,因此對AIMD平衡軌跡中的構(gòu)型進行篩選并剔除重復構(gòu)型,最終共計選出871 個結(jié)構(gòu)用于構(gòu)建深度學習勢函數(shù),隨機選取其中88 個(占比10%)結(jié)構(gòu)作為驗證集,剩余783 個(占比90%)結(jié)構(gòu)用于勢函數(shù)的訓練.

圖2 β-Si3N4 的晶體結(jié)構(gòu) (a) 用于第一性原理分子動力學計算的體系(原子個數(shù)為140 個);(b) 用于深度學習分子動力學模擬的體系(原子個數(shù)為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 深度學習勢函數(shù)訓練

本文采用DeepMD-kit[21]軟件包訓練深度學習勢,該軟件基于當下流行的深度學習框架Tensorflow[22]進行二次開發(fā),能夠?qū)崿F(xiàn)多體勢訓練并保持計算精度幾乎不損失.依據(jù)原子神經(jīng)網(wǎng)絡勢的理論框架[12,23]可知,對于包含N個原子的體系,體系總能量可歸結(jié)為每個原子的能量之和:

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

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

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

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

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

3 結(jié)果與討論

3.1 勢函數(shù)訓練及結(jié)果驗證

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

圖3 (a) 第一性原理計算體系能量和深度學習勢計算體系能量對比關(guān)系圖;(b)—(d) 第一性原理計算原子受力和深度學習勢計算原子受力對比關(guān)系圖.其中圖中直線代表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 方法計算了體系的徑向分布函數(shù)進行對比,來驗證深度學習勢在描述結(jié)構(gòu)演化過程中的精度和準確性.在相同溫度條件下以AIMD 初始結(jié)構(gòu)進行深度學習分子動力學模擬,計算體系原子數(shù)目均為140.兩種模擬得到的Si-Si,Si-N,N-N 的徑向分布函數(shù)如圖4 所示,可以明顯看到深度學習勢準確地復現(xiàn)了AIMD 徑向分布函數(shù)峰的形狀和位置,表明深度學習勢函數(shù)可以很好地描述體系的動力學性質(zhì).

圖4 (a)—(c) 300 K 溫度下DPMD 與AIMD 計算徑向分布函數(shù)對比圖;(d)—(f) 1000 K 溫度下DPMD 與AIMD 模擬徑向分布函數(shù)對比圖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 弛豫時間及微波介電函數(shù)計算

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

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

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

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

其中t為時間.

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

圖5 運用深度學習分子動力學計算的700 K 溫度下氮化硅偶極矩自相關(guān)函數(shù) (a) 偶極矩自相關(guān)函數(shù)隨時間的變化;(b) 偶極矩自相關(guān)極值點對數(shù)值隨時間的變化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.

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

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

圖6 (a) β-Si3N4 在不同溫度時的頻率相關(guān)介電函數(shù)實部;(b) β-Si3N4 在不同溫度時的頻率相關(guān)介電函數(shù)虛部.其中文獻[34]僅給出了實部值,未給出虛部值.文獻[35]給出了8—12 GHz 頻率范圍內(nèi)的介電函數(shù)實驗測量值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 勢來說,其構(gòu)造過程是先選擇合適的數(shù)學表達式,然后通過實驗數(shù)據(jù)擬合確定函數(shù)中的參數(shù),因此經(jīng)驗勢只能描述已知的性質(zhì).因此,利用Tersoff 勢來計算高溫氮化硅的介電特性還有待商榷.相比之下,基于第一性原理計算求解電子結(jié)構(gòu)所得到結(jié)果是值得信賴的,而深度學習勢是基于第一性原理數(shù)據(jù)進行訓練擬合,因此深度學習勢計算得到的氮化硅介電函數(shù)結(jié)果更加合理.

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

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

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

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

圖8 DPMD (Nvidia RTX 3080 GPU 計算)與AIMD(48 個Intel Xeon Platinum 9242 CPU 計算)計算速度結(jié)果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 結(jié)論

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