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

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

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

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

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

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

1 固體界面摩擦性能的計(jì)算方案

1.1 自動(dòng)化計(jì)算流程

固體界面摩擦性能整體設(shè)計(jì)流程如圖1所示，具體如下：(1)讀取界面結(jié)構(gòu)文件；(2)分析界面對(duì)稱性，在x、y和z方向上移動(dòng)界面構(gòu)建計(jì)算摩擦性能所需的所有位置文件，命名為POSCAR-i-j-k(i、j和k是界面在x、y和z方向上移動(dòng)的步數(shù)).根據(jù)界面結(jié)構(gòu)文件構(gòu)建POTCAR、KPOINTS以及INCAR文件，新建文件夾i-j-k，將對(duì)應(yīng)的POSCAR-i-j-k以及VASP計(jì)算所需的其他文件復(fù)制到該文件夾；(3)依次進(jìn)入i-j-k文件夾并調(diào)用VASP執(zhí)行總能計(jì)算；(4)依次進(jìn)入i-j-k文件夾，提取所有計(jì)算構(gòu)型的總能，存入Energy[i][j][k]三維數(shù)組.對(duì)每組(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é)性能計(jì)算流程

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

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

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

1.3 第一性原理計(jì)算方法

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

1.4 摩擦性能的計(jì)算方法

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

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

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

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

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

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

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

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

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

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

2 結(jié)果與討論

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

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

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

Fig.4 The system energy as a function of interplanar distance z for high-symmetry sites in graphene/graphene system圖4 石墨烯/石墨烯滑動(dòng)體系高對(duì)稱位置的能量與界面距離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)系擬合精度，對(duì)z方向30個(gè)數(shù)據(jù)點(diǎn)進(jìn)行分組，每組包含8個(gè)數(shù)據(jù)點(diǎn)，最后1組數(shù)據(jù)包含6個(gè)數(shù)據(jù)點(diǎn)，然后采用最小二乘法對(duì)每組數(shù)據(jù)進(jìn)行擬合.E(z)的函數(shù)形式采用12階多項(xiàng)式，如式(1)所示.最后，檢查所有擬合數(shù)據(jù)，結(jié)果表明，采用這種方法擬合的數(shù)據(jù)絕對(duì)誤差小于1.0×10-5eV.圖5(a)給出了其中1組數(shù)據(jù)和擬合曲線，圖5(b)給出了數(shù)據(jù)擬合誤差.

2.3 勢(shì)能面PES

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

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

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

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

2.5 該方法的局限性

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

3 結(jié)論

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

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

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

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