基于第一性原理的固體界面摩擦學(xué)性能高通量計算研究

何文豪, 魯志斌*, 范曉麗, 周 峰, 劉維民

(1.中國科學(xué)院蘭州化學(xué)物理研究所 固體潤滑國家重點(diǎn)實(shí)驗(yàn)室, 甘肅 蘭州 730000;2.西北工業(yè)大學(xué) 材料學(xué)院 先進(jìn)潤滑與密封材料中心 凝固技術(shù)國家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710000)

摩擦學(xué)綜合了物理學(xué)、化學(xué)、材料科學(xué)和力學(xué)等學(xué)科領(lǐng)域，主要研究相對運(yùn)動或有相對運(yùn)動趨勢的界面間的摩擦、潤滑和磨損行為及其機(jī)理.影響界面摩擦性能的因素是多方面的，主要包括界面性質(zhì)、負(fù)載、服役環(huán)境和服役溫度等.正是因?yàn)槿藗冃枰C合考慮各種因素來理解摩擦，從非常小尺度下的量子相互作用(界面的電子結(jié)構(gòu)、吸附等)到大尺度下的宏觀變量(負(fù)載、服役溫度等)，這使得摩擦學(xué)的研究幾乎不可能沒有計算機(jī)模擬技術(shù)的輔助.在過去的幾十年里，隨著計算機(jī)硬件的持續(xù)發(fā)展，計算機(jī)模擬技術(shù)逐漸成為推動摩擦學(xué)發(fā)展的關(guān)鍵技術(shù)之一.

計算機(jī)模擬方法主要包括基于密度泛函理論(DFT)的第一性原理方法、分子動力學(xué)方法(MD)和有限元分析(FEA)方法等.這些方法在不同的尺度上對理解材料性能都起到了重要作用.盡管對摩擦性能的研究來說，DFT方法限制了界面結(jié)構(gòu)(局限于無限平坦的界面)，但它可以深入到量子尺度，從界面電子結(jié)構(gòu)的角度理解摩擦[1-5].據(jù)此，作者所在課題組提出了一種實(shí)現(xiàn)超滑的新策略-壓力誘導(dǎo)超滑[6].科研人員也通過DFT方法解釋了許多摩擦試驗(yàn)現(xiàn)象，比如界面鈍化對金剛石[7-8]、類金剛石[9]、MoS2[10]和 Fe[11-13]摩擦 性 能的影響.這使得DFT方法成為理解摩擦性能和解釋摩擦試驗(yàn)現(xiàn)象的重要技術(shù)手段.

近年來，在材料科學(xué)領(lǐng)域，高通量DFT計算模擬研究越來越普遍[14].蘭州化學(xué)物理研究所周等[15]提出了離子液體和酯類化合物兩類潤滑劑物性參數(shù)和摩擦系數(shù)的高通量分子動力學(xué)計算框架.意大利摩德納大學(xué)的Pizzi等[16]基于“自動化交互基礎(chǔ)設(shè)施和數(shù)據(jù)庫的計算科學(xué)框架”(AiiDA)，提出了第一個固體界面摩擦學(xué)性能的高通量DFT計算方法[17].盡管這個方法實(shí)現(xiàn)了固體界面的勢能面(PES)和理想剪切強(qiáng)度(τ)的計算，但是該方法僅適用于零負(fù)載下PES和τ的計算，不適用于任意負(fù)載(Fz).而寬Fz范圍內(nèi)摩擦性能的計算更能體現(xiàn)實(shí)際情況，更有利于解釋試驗(yàn)現(xiàn)象.另外，通過建立量子尺度摩擦性能與負(fù)載間的關(guān)系，更有助于理解量子尺度摩擦的物理本質(zhì).

目前，特定負(fù)載下固體界面摩擦性能的DFT計算尚未實(shí)現(xiàn)自動化和高通量，計算模型構(gòu)建以及數(shù)據(jù)后處理浪費(fèi)了科研人員的大量時間.通常要完整計算1個固體界面體系的摩擦性能，首先需要建立約500～1 000個計算模型，然后準(zhǔn)備DFT計算軟件所需的輸入文件，提交計算任務(wù)，計算結(jié)束后手動搜索每個構(gòu)型的總能，擬合能量與界面距離的關(guān)系，然后才能獲得PES，最后在PES上尋找能量最低的滑動路徑，計算摩擦力以及摩擦系數(shù).整個過程人工操作耗費(fèi)的時間約60小時.另外，科研人員多研究1個體系，耗費(fèi)的時間就增加1倍.因此，非常有必要提出1個能夠?qū)崿F(xiàn)任意負(fù)載固體界面摩擦學(xué)性能高通量DFT計算方法.

鑒于此，本文作者提出了1個能夠?qū)崿F(xiàn)自動化建模、計算任務(wù)自動提交管理和智能化數(shù)據(jù)后處理的方法，即自動提取計算結(jié)果、擬合數(shù)據(jù)、繪制勢能面、搜索最優(yōu)滑動路徑、輸出摩擦力和摩擦系數(shù)的固體界面摩擦性能高通量DFT計算方法，該方法能夠極大地節(jié)約科研人員研究固體界面摩擦性能所需的時間.

1 固體界面摩擦性能的計算方案

1.1 自動化計算流程

固體界面摩擦性能整體設(shè)計流程如圖1所示，具體如下：(1)讀取界面結(jié)構(gòu)文件；(2)分析界面對稱性，在x、y和z方向上移動界面構(gòu)建計算摩擦性能所需的所有位置文件，命名為POSCAR-i-j-k(i、j和k是界面在x、y和z方向上移動的步數(shù)).根據(jù)界面結(jié)構(gòu)文件構(gòu)建POTCAR、KPOINTS以及INCAR文件，新建文件夾i-j-k，將對應(yīng)的POSCAR-i-j-k以及VASP計算所需的其他文件復(fù)制到該文件夾；(3)依次進(jìn)入i-j-k文件夾并調(diào)用VASP執(zhí)行總能計算；(4)依次進(jìn)入i-j-k文件夾，提取所有計算構(gòu)型的總能，存入Energy[i][j][k]三維數(shù)組.對每組(i,j)，擬合總能與界面距離z的關(guān)系E(z).采用最小二乘法擬合總能與界面距離關(guān)系E(z).E(z)的函數(shù)形式采用12階多項(xiàng)式，如式(1)所示：

Fig.1 The operation process of the solid interface tribological properties圖1 固體界面摩擦學(xué)性能計算流程

(5)求解不同負(fù)載Fz對 應(yīng)的能量，繪制Fz對應(yīng)的PES；(6)搜索PES中能壘最低的滑動路徑，計算相應(yīng)的摩擦力f和摩擦系數(shù)μ.

1.2 固體界面摩擦性能高通量設(shè)計方案

對于固體界面摩擦性能第一性原理計算，計算量的大小取決于界面沿著x，y和z方向平移后，產(chǎn)生的計算模型數(shù)量，通常來說這個數(shù)量在1 000左右，逐個計算這些模型的能量非常耗時，因此，有必要通過高通量的方法并行-并發(fā)計算這些構(gòu)型的能量.本文中使用圖2所示的設(shè)計思路實(shí)現(xiàn)固體界面摩擦性能的高通量計算.首先基于界面結(jié)構(gòu)生成計算摩擦性能所需的約1 000個計算模型，然后通過調(diào)用服務(wù)器的作業(yè)調(diào)度系統(tǒng)，同時將這些計算任務(wù)提交到服務(wù)器，實(shí)現(xiàn)計算任務(wù)的并發(fā)處理，并調(diào)用Vienna Ab initio Simulation Package(VASP)軟件，對同一任務(wù)使用多核并行計算，實(shí)現(xiàn)1個固體界面體系摩擦性能的并發(fā)并行計算.

1.3 第一性原理計算方法

本文中第一性原理的總能計算利用基于DFT[18]的Vienna Ab-initio Simulation Package (VASP)軟件包[19-21].電子與離子實(shí)之間的相互作用采用贗勢方法處理.贗勢使用投影綴加平面波(PAW)方法[22]，電子之間的交互關(guān)聯(lián)勢使用廣義梯度近似(GGA)中的Perdew-Burke-Ernzerhof (PBE)[23]泛函.該高通量計算方法中，范德華相互作用被設(shè)置為可選變量，需要在程序輸入文件中指定.對于文中石墨烯/石墨烯滑動界面的計算，修正格林函數(shù)用于修正范德華力(vdW)色散相互作用[24].該高通量計算方法會通過總能收斂性測試確定平面波展開的截斷能量，收斂性測試表中是能量差小于1 meV/atom.電子步的能量收斂標(biāo)準(zhǔn)是1×10-6eV.布里淵區(qū)使用Monknorst-Pack的K點(diǎn)網(wǎng)格(k-mesh)，程序會通過總能收斂性測試確定K點(diǎn)大小.在總能計算過程中，不弛豫原子的位置，所有體系都使用周期性邊界條件.

1.4 摩擦性能的計算方法

首先計算固體界面滑動系統(tǒng)的勢能面PES，Zhong等[25]提出的固體界面摩擦過程中的勢能V定義如下：

Fig.2 The high throughput algorithm of the solid interface tribological properties圖2 固體界面摩擦學(xué)性能高通量設(shè)計方案

其中Ead是界面能，x和y是界面在這兩個方向上移動的距離，z是界面距離，F(xiàn)z是垂直界面方向的負(fù)載，V0是PES中最低的勢能，定義為

其中x0和y0是能量最低的界面在x和y方向上移動的距離.界面能Ead的定義如下：

其中EInterface是界面體系的能量，Eup和Edown分別是界面上部分和下部分的能量.負(fù)載Fz與吸附能Ead之間的關(guān)系滿足下面的方程：

其中A是界面面積，本文中A定義為界面的表觀接觸面積.摩擦力f定義如下：

其中 ΔV是PES勢能面滑動路徑上鞍點(diǎn)能量與最低點(diǎn)能量之間的差值，d是鞍點(diǎn)與最低點(diǎn)的距離.根據(jù)方程(7)，本文中將摩擦力轉(zhuǎn)化為界面剪切強(qiáng)度τ：

摩擦系數(shù) μ定義為

1.5 等負(fù)載摩擦計算的實(shí)現(xiàn)方法

盡管VASP軟件可以計算不同應(yīng)力下體系總能，但是對于界面體系摩擦性能的計算，這種方法存在兩方面的局限.其一，摩擦性能計算需要僅在垂直于界面的方向施加載荷，但是VASP只能在3個方向上同時施加相同的載荷；其二，考慮到周期性邊界條件的使用，界面體系需要包含真空層，而VASP軟件無法直接對包含真空層的體系施加載荷.因此，該高通量計算方法沒有使用VASP直接施加載荷，而是通過Zhong等[23]給出的數(shù)值擬合的方法計算等載荷下的摩擦性能.具體實(shí)現(xiàn)流程如下：(a)計算不同滑移路徑和界面距離對應(yīng)的體系總能.然后對每條滑移路徑上的每個位移點(diǎn)，執(zhí)行b～e步驟；(b)通過方程(1)擬合界面總能與界面距離之間的關(guān)系；(c)使用方程(5)計算負(fù)載與界面距離之間的定量關(guān)系；(d)通過方程(5)計算給定負(fù)載對應(yīng)的界面距離，根據(jù)這個界面距離，使用方程(1)計算給定負(fù)載對應(yīng)的界面總能；(e)通過方程(2)計算給定負(fù)載體系的勢能；(f)根據(jù)給定負(fù)載對應(yīng)的滑動過程中所有體系的勢能繪制等負(fù)載的PES，根據(jù)這個PES計算給定負(fù)載下的界面剪切強(qiáng)度和摩擦系數(shù).

2 結(jié)果與討論

本文中以石墨烯/石墨烯滑動體系為例，測試了提出的固體界面摩擦學(xué)性能高通量第一性原理計算方法的可靠性.界面在x、y和z方向上分別移動的步數(shù)為12、12和30，考慮到界面結(jié)構(gòu)的對稱性，總共構(gòu)建了1 080個計算模型.x和y方向每步移動的距離約為0.21 ?，在z方向每步移動的距離為0.10 ?.這些模型中包含石墨烯/石墨烯滑移體系全部高對稱位置，包括TOP(T)位置、Bridge(B)位置、Hollow(H)位置以及Saddle(S)位置，如圖3所示.在1 024核服務(wù)器上運(yùn)行該計算軟件，每個計算任務(wù)使用4個核，服務(wù)器同時并行運(yùn)算256個計算任務(wù).

2.1 能量與界面距離關(guān)系

計算完成后，提取了所有體系的總能.圖4給出了石墨烯/石墨烯滑動體系高對稱位置的能量與界面距離z之間的關(guān)系.顯然，隨著z增大，體系總能降低，當(dāng)z約為界面間的平衡距離時，界面能降到最低，之后體系總能隨著z增大而增大.值得注意的是，總能最低時，不同位置對應(yīng)的界面距離存在小差異，Top位置對應(yīng)的z約 為3.45 ?，其他3個位置對應(yīng)的z約為3.30 ?，如圖4所示.另外，當(dāng)界面距離z約為1.63 ?時，4個高對稱位置對應(yīng)的體系總能近似相等，如圖4所示.說明z為1.63 ?時，界面滑動需克服的能壘約為零，從而可能實(shí)現(xiàn)超滑.

Fig.3 The high-symmetry Top (T), Bridge (B), Hollow (H) and Saddle (S) sites in the Graphene/Graphene sliding system圖3 石墨烯/石墨烯滑動體系的高對稱位置

Fig.4 The system energy as a function of interplanar distance z for high-symmetry sites in graphene/graphene system圖4 石墨烯/石墨烯滑動體系高對稱位置的能量與界面距離z之間的關(guān)系

2.2 能量與界面距離關(guān)系擬合

之后，擬合了體系能量E與界面距離z之間的關(guān)系E(z).E(z)關(guān)系的擬合精度決定了給定負(fù)載Fz求解能量的準(zhǔn)確度，從而決定了摩擦力f和摩擦系數(shù)μ的誤差.因此，E(z)關(guān)系需要非常高的擬合精度，為了提高E(z)關(guān)系擬合精度，對z方向30個數(shù)據(jù)點(diǎn)進(jìn)行分組，每組包含8個數(shù)據(jù)點(diǎn)，最后1組數(shù)據(jù)包含6個數(shù)據(jù)點(diǎn)，然后采用最小二乘法對每組數(shù)據(jù)進(jìn)行擬合.E(z)的函數(shù)形式采用12階多項(xiàng)式，如式(1)所示.最后，檢查所有擬合數(shù)據(jù)，結(jié)果表明，采用這種方法擬合的數(shù)據(jù)絕對誤差小于1.0×10-5eV.圖5(a)給出了其中1組數(shù)據(jù)和擬合曲線，圖5(b)給出了數(shù)據(jù)擬合誤差.

2.3 勢能面PES

圖6示出了不同負(fù)載對應(yīng)的勢能面(PES)文件，如圖6(a)所示，0 GPa條件下，PES比較平坦，T位置勢能最高，H位置勢能最低，界面沿著H-S-H滑動.隨著負(fù)載增大，PES起伏逐漸變大，界面滑動需克服的能壘也逐漸增大.當(dāng)負(fù)載約為180 GPa時，PES起伏最大，滑動能壘也達(dá)到最大，如圖6(c)所示.之后隨著負(fù)載進(jìn)一步增大，PES起伏開始變小，滑動能壘也開始變小，當(dāng)負(fù)載約為253 GPa時，PES變得平坦，此時界面滑移能壘與0 GPa時近似相等.隨著負(fù)載進(jìn)一步增大，PES起伏又逐漸變大，但此時能量最低的位置變?yōu)門位置，H位置轉(zhuǎn)變?yōu)槟芰孔罡叩狞c(diǎn)，即出現(xiàn)了PES反轉(zhuǎn).本文中得出的PES隨著負(fù)載的變化趨勢與課題組通過傳統(tǒng)的方法計算的結(jié)果一致[6]，這驗(yàn)證了該方法的可靠性.

2.4 界面剪切強(qiáng)度和摩擦系數(shù)

圖7示出了計算的石墨烯/石墨烯滑動體系的界面剪切強(qiáng)度 τ、摩擦系數(shù) μ與負(fù)載Fz之間的關(guān)系.可以看出，隨著Fz的 增大，τ先增大然后減小，最后增大.當(dāng)Fz為 0 GPa時，τ最小，約0.1 GPa.當(dāng)Fz約為180 GPa時，τ達(dá)到極大值，約3.9 GPa.當(dāng)Fz約 為253 GPa時，τ降到極小值，約為0.3 GPa，出現(xiàn)了壓力誘導(dǎo)超滑現(xiàn)象.這個結(jié)果與PES隨著負(fù)載的變化結(jié)果一致.也與課題組通過傳統(tǒng)的方法計算的定性結(jié)果[6]一致，但在定量比較時出現(xiàn)了偏差.傳統(tǒng)方法給出負(fù)載約300 GPa時出現(xiàn)了壓力誘導(dǎo)超滑現(xiàn)象[6]，但本文中計算的壓力誘導(dǎo)超滑的負(fù)載值約為253 GPa，這主要是因?yàn)橛嬎惴椒ǖ牟煌?文獻(xiàn)[6]中使用了等界面距離的計算方法，而本文中使用了等負(fù)載的計算方法.本文中使用的等負(fù)載方法的計算結(jié)果更有利于與試驗(yàn)進(jìn)行直接對比，摩擦系數(shù) μ隨著負(fù)載的變化曲線也觀察到了壓力誘導(dǎo)超滑現(xiàn)象，如圖7(b)所示.當(dāng)負(fù)載為50～200 GPa時， μ約為0.02～0.03，這與試驗(yàn)值約1 0-2一致[26].

Fig.5 (a) The relationship between energies of system and interface-distances (the black blocks are the calculating energies and the blue line is the fitting curve) and (b) errors between calculating energies and fitting energies圖5 (a)體系能量E與界面距離z的關(guān)系(黑色方塊表示計算數(shù)據(jù)，藍(lán)色實(shí)線表示擬合曲線)；(b)計算和擬合數(shù)據(jù)之間的誤差

Fig.6 The potential energy surface corresponding to different loads in the graphene/graphene sliding system圖6 石墨烯/石墨烯滑動體系中不同負(fù)載對應(yīng)的勢能面

2.5 該方法的局限性

盡管該高通量計算方法實(shí)現(xiàn)了固體界面摩擦學(xué)性能計算的自動化以及高通量計算，并使用石墨烯/石墨烯滑動體系驗(yàn)證了該方法的可靠性，但該方法仍然存在一些局限性.首先，由于DFT計算本身無法直接引入溫度，因此該方法也沒有考慮溫度對摩擦性能的影響；其次，該方法不僅能夠計算范德華相互作用結(jié)合的簡單二維材料，也可以計算氧化物和金屬之間的同/異質(zhì)界面，但需要注意的是，該方法沒有考慮摩擦過程中界面之間原子的交換以及界面摩擦化學(xué)反應(yīng).摩擦過程中非晶以及較粗糙的界面可能更容易發(fā)生界面原子的轉(zhuǎn)移，因此，使用該方法計算非晶或者較粗糙的界面時，需要仔細(xì)考慮計算模型以及結(jié)果的合理性；另外，該方法沒有考慮氣氛環(huán)境的影響，氣氛環(huán)境對固體界面摩擦的影響普遍存在，因此，非常有必要在該方法中引入氣氛環(huán)境的影響，這也是該方法下一步需要改進(jìn)的方向.

Fig.7 (a) Interfacial shear strength τ and (b) friction coefficient μ as a function of the load Fz in the graphene/graphene sliding system圖7 石墨烯/石墨烯滑動體系中(a)界面剪切強(qiáng)度τ和(b)摩擦系數(shù)μ與負(fù)載Fz之間的的關(guān)系

3 結(jié)論

報道了1個基于第一性原理的固體界面摩擦學(xué)性能高通量計算方法，并以石墨烯/石墨烯滑動體系為例驗(yàn)證了該方法的可靠性.主要結(jié)論如下：

a.該方法實(shí)現(xiàn)了固體界面摩擦學(xué)性能計算的自動化，即自動構(gòu)建計算模型、自動提交管理計算任務(wù)和智能分析計算數(shù)據(jù)，直接輸出界面的勢能面、摩擦力和摩擦系數(shù)，能夠極大地節(jié)約科研人員使用第一性原理研究固體界面摩擦性能所需的時間.

b.該方法實(shí)現(xiàn)了固體界面摩擦性能計算的高通量，即1 000多個計算任務(wù)的并發(fā)并行計算，提高了固體界面摩擦性能的計算效率.

c.使用該方法計算的石墨烯/石墨烯滑動體系的勢能面與使用傳統(tǒng)方法的計算結(jié)果一致，摩擦系數(shù)與試驗(yàn)結(jié)果符合，從而驗(yàn)證了該方法的可靠性.