利用雙晶格勢進行硅的分子動力學模擬及其有效性

魏崇陽,張 輝,劉仲武,鐘喜春,焦東玲,邱萬奇,余紅雅

(華南理工大學 材料科學與工程學院,廣東 廣州 510640)

1 前 言

硅(Si)是重要的半導體元素[1-8]。在過去的幾十年中,基于實驗理論的模擬計算方法在Si及其化合物的研究中發(fā)揮了重要作用[9]。隨著研究對象尺度的不斷減小,微/納米尺度Si團簇的物理化學性質的實驗研究變得越來越困難,原子尺度的計算模擬變得尤為重要。由于Si是共價晶體,原子之間由共價鍵結合,在原子尺度的計算中對共價鍵的描述十分困難。從理論上講,密度泛函理論(DFT)計算可以給出最精確的計算結果,但即使是由數(shù)百個原子組成的系統(tǒng)也需要大量的計算時間以及非常高的計算能力,這使得量化計算適用范圍受限。因此,在研究稍大系統(tǒng)的性質時,最常用的是基于相互作用經驗勢的分子動力學(MD)模擬方法,通過構造經驗勢,并且應用于分子動力學模擬中來計算系統(tǒng)的物性參數(shù)。

目前,Si的原子間作用勢包括Tersoff勢以及其修正勢[10-13]、Stillinger-Weber(SW)勢[14]、依賴于環(huán)境的原子間作用(EDI)勢[15]、電荷優(yōu)化多體(COMB)勢[16]、修飾嵌入原子(MEAM)勢[17]及其修正勢等[18]。其中,Tersoff勢、SW 勢、EDI勢和COMB 勢是多體勢。在這些作用勢中,除了描述兩個原子之間的相互作用的對勢之外,還存在一個三體勢,用以描述三個原子之間的相互作用或定義所謂的鍵角。此外,COMB勢還考慮了電荷的影響。這些多體勢可以為Si的某些物理性質提供適當?shù)拿枋?。許多研究者[10-18]通過計算基于不同作用勢的Si的相關參數(shù),例如彈性常數(shù),結合能和徑向分布函數(shù)等,來驗證作用勢的有效性。雖然這些物理參數(shù)反映了物質本身的某些屬性,但最直接的驗證標準應該是應用這些原子間作用勢能否在分子動力學模擬中獲得Si的金剛石結構。

研究表明,即使應用一個與實驗數(shù)據(jù)非常符合的原子間作用勢,在分子動力學模擬中也無法形成金剛石結構。例如,應用描述金屬元素和合金中的嵌入原子(EAM)勢[19],在分子動力學模擬中可以得到與實驗一致的金屬晶格結構[20-21]。然而在分子動力學模擬中應用上述原子間作用勢,模擬系統(tǒng)卻不能得到Si的金剛石結構,即上述作用勢無法通過Si的晶體結構的驗證測試。因此,本研究為Si構建了一個新的原子間作用勢,利用該作用勢在分子動力學模擬中得到了Si的金剛石結構,并據(jù)此獲得了相關的物理參數(shù),用以驗證該作用勢的有效性。

金剛石結構可以看作是兩個面心立方(fcc)晶格的組合,其晶胞中有八個原子,因此可以從兩個fcc晶格原子間作用勢和兩個晶格之間的原子間作用勢出發(fā)構建金剛石結構。目前,用于fcc晶格的模擬主要為EAM 勢[19-21],但是在EAM 勢中引入的參數(shù)過多,很難定義兩個fcc晶格之間的原子間作用。這里,僅考慮簡單的Lenard-Jones(LJ)勢。研究表明,在分子動力學模擬中LJ勢的基態(tài)對應于六方密排(hcp)晶格(c/a≈1.633,a和c為hcp結構的晶格常數(shù)),這與計算結果一致[22-23]。在某些模擬情況下,系統(tǒng)還可以顯示亞穩(wěn)態(tài)fcc晶格[24]。將模擬結果與理想的hcp 和fcc晶格比較,可以鑒別出模擬系統(tǒng)的布拉格點陣。進一步的研究表明,在模擬系統(tǒng)中不設置任何初始布拉格點陣,利用幾種LJ勢的組合,構造出的雙晶格(DL)勢可以得到金剛石和石墨結構[25]以及更復雜的鈣鈦礦ABO3結構[26]。這些結果有助于理解和重新思考LJ勢的重要性以及晶體結構的形成。

本研究在分子動力學模擬中通過構造雙晶格勢來獲得Si的金剛石結構,并應用雙晶格勢獲得了Si的結晶溫度、晶格常數(shù)和彈性常數(shù)等物理參數(shù),并與其他原子間相互作用勢和實驗測試結果進行了比較。

2 模型與模擬

2.1 模型

在雙晶格作用勢中,通過將兩個fcc晶格原子間勢和兩個晶格之間的原子間勢組合來描述晶體結構。用于描述fcc晶格的LJ勢可表示為:

式中:ε是勢阱的深度,σ原子間相互作用勢為零的有限距離,r是原子之間的距離。

首先建立一個多原子相互作用的模型,在該模型中引入兩種類型的原子,分別表示為S1和S2,它們在物理屬性上彼此相同,對應于Si金剛石結構原子的兩個自旋結構。例如,S1類型原子對應四個自旋向上,S2類型原子對應四個自旋向下。S1和S2原子的晶格和晶格常數(shù)相同。在其他模型中,所有Si原子都相同,而該模型中的S1和S2原子類型被視為兩個可區(qū)分的Si原子類型。S1和S2晶格之間存在相互作用。在這里,LJ勢應用于S1晶格,S2晶格以及S1和S2晶格之間的相互作用,原子之間的平衡距離和結晶溫度由LJ勢參數(shù)確定。

2.2 模擬過程

模擬過程中,在不設置任何初始布拉格點陣情況下,在模擬盒子中平均隨機創(chuàng)建S1和S2原子。應用周期性邊界條件,對于LJ勢,截止距離表示為rc。系統(tǒng)的最高溫度不高于溫度T0(4 000 K),然后執(zhí)行能量最小化過程。在T0溫度下執(zhí)行NPT 系綜操作,時長為1 000 ps,時間步長為10-3ps。系統(tǒng)的溫度從T0開始降低,步長為n。在每個溫度T處執(zhí)行一次NPT系綜操作,持續(xù)時間為100～1 000 ps。在整個模擬過程中壓力始終保持為零。表1列出了分子動力學模擬中LJ勢的參數(shù)。

表1 分子動力學模擬中LJ勢的參數(shù)Table 1 Parameters for LJ potentials in MD simulations

為了獲得基于DL勢的Si晶體應力與應變關系,在模擬系統(tǒng)中創(chuàng)建理想的Si晶體結構。S1和S2原子的初始fcc晶格具有相同的晶格常數(shù)(a=5.341?),并且S1晶格相對于S2晶格在x,y和z方向上的相對位移分別為0.25、0.25 和0.25 個單位長度。在300 K的溫度下執(zhí)行NPT 系綜操作,時長為1 000 ps。使晶體發(fā)生微小變形,之后計算由變形引起的應力值,并獲得應力與應變的關系。由此,可以計算出彈性常數(shù)c11、c12和c44的值,計算代碼參考文獻[27]。

本研究還將基于DL勢計算獲得的物理參數(shù)與其他原 子 相 互 作 用 勢(Tersoff、SW、EDI、COMB 和MEAM 勢)的結果進行了比較。本次模擬的系統(tǒng)原子數(shù)分別為512、1 000和8 000。分子動力學模擬使用的軟件為lammps[28],并通過VESTA[29]軟件進行可視化分析。

2.3 模擬系統(tǒng)晶體結構識別

檢查識別模擬系統(tǒng)顯示為Si的金剛石結構是測試DL勢有效性的首要條件。但是lammps軟件缺乏識別模擬系統(tǒng)晶體結構的功能,因此本研究計算了模擬系統(tǒng)中原子之間的距離分布函數(shù)和原子與其最近鄰原子之間的角度的分布函數(shù),即ρ(d)和ρ(θ)。從lammps軟件獲得理想Si晶體結構,其分布函數(shù)由少量離散值組成。如果模擬的分布函數(shù)在這些相應位置顯示為非零值,則可認為該系統(tǒng)具有與理想Si晶體相同的晶體結構,從而檢查模擬系統(tǒng)與參考系統(tǒng)之間的差異以進行晶格識別。在模擬中,得到任一原子在特定時間和溫度下的(x,y,z)坐標,可以計算出原子之間的距離及原子與其最鄰近原子之間的角度。分布函數(shù)ρ(d)(或ρ(θ))表示距離d(或角度θ)在d-d+dd(或θ-θ+dθ)范圍內的計數(shù)值或強度。在計算中,如果一個原子R及其最近鄰原子之間的距離為dR,那么一個原子的最近鄰原子數(shù)始終大于1,這些距離dR的最小距離為dm,并且dR的值不同。如果dR/dm＜1.1,則該原子可被視為原子R的最近鄰原子之一。

此外,分布函數(shù)ρ(d)與徑向分布函數(shù)g(r)相似,但是它們的計算方法不同。角度θ是某些晶體結構(如金剛石結構)的鍵角。對于Si的晶格識別,為了計算分布函數(shù),可將S1和S2原子視為同一種原子,也可將其單獨用作S1(或S2)子系統(tǒng)。

3 結果與分析

3.1 晶體相變和結構識別

本研究根據(jù)DL 勢和其他原子間勢(Tersoff、MEAM、EDI、SW 和COMB勢)計算了不同溫度下Si原子系統(tǒng)的能量和體積。某些勢會有一些修正版本,例如Tersoff的Tersoff-Mod(TM)勢[12]和Tersoff-Mod C(TMC)勢[13],以及MEAM 的MEAM-Spline(MEAMS)勢[18]。圖1 顯示了每個原子的能量和每個原子的體積與溫度的關系。圖1(a)表明,對于所有作用勢,系統(tǒng)的能量隨著溫度的降低而降低。特別是對于Tersoff勢,模擬系統(tǒng)始終處于簡單的結晶狀態(tài),一旦溫度超過1 800 K,系統(tǒng)的體積就會迅速增加到無窮大,導致模擬失敗。這意味著具有Tersoff勢的系統(tǒng)將永遠不會進入液態(tài),并且始終顯示金剛石結構。對于TM 勢,系統(tǒng)可以進入液態(tài)。隨著溫度的降低,每個原子的能量不斷減少,直到溫度接近1 300 K,并且能量突然改變,表明存在相變。但是,在使用EDI、MEAM、MEAMS、SW 和COMB勢進行的模擬中,系統(tǒng)體積隨溫度降低突然增加(見圖1(b))。對于基于DL勢的模擬還存在一個明顯的熱滯。從高溫冷卻時,能量和體積都會降低,并且存在液態(tài)晶態(tài)相變,這與實驗現(xiàn)象一致。在DL 勢的情況下,每個原子的能量在100 K 下約為2.18 eV,僅為其他作用勢的一半。但是,該能量仍可以確保結晶溫度與實驗結果符合,體積值也與實驗數(shù)據(jù)非常吻合(見表2)。

圖1 基于不同作用勢單原子的能量(a)和體積(b)與溫度的關系Fig.1 Dependence of energy per atom(a)and volume per atom(b)on the temperature with different interatomic potentials applied

表2 基于不同原子間作用勢獲得的結晶溫度T c、晶格常數(shù)a 和彈性常數(shù)(c 11,c 12,c 44)Table 2 Crystallization temperature T c,lattice constant a,and elastic constants(c 11,c 12,and c 44)obtained with different interatomic potentials applied

圖2和3分別顯示在T=100 K 和4 000 K 時,基于不同原子間作用勢模擬得到的原子之間的距離ρ(d)和原子與其最近鄰原子之間的角度ρ(θ)的分布函數(shù)。圖2中,對于某些確定的d和θ,理想的金剛石結構的ρ(d)和ρ(θ)不為零。當系統(tǒng)冷卻后,Tersoff勢的ρ(d)和ρ(θ)與理想Si晶體結構的結構相匹配,表明Tersoff勢有利于保持金剛石結構的穩(wěn)定,而具有TM 勢的系統(tǒng)顯示出明顯的不同。例如,在T=100 K時,ρ(θ)在109.5°處雖然有一個弱峰,但如圖2(a)所示,它卻不是晶態(tài)的,表明此刻Tersoff勢模擬系統(tǒng)處于結晶狀態(tài)和液態(tài)的混合狀態(tài)。在T=4 000 K 時,所有的峰消失,所有系統(tǒng)都處于液態(tài)(見圖3(a))。

圖2 基于不同作用勢原子之間的距離ρ(d)和原子與其最近鄰原子之間的角度ρ(θ)的分布函數(shù)(T=100 K;COMB勢T=50 K)Fig.2 Distribution functions of the distances between atoms and the angles between the lines linking one atom with its nearest neighborsρ(d) (a)andρ(θ) (b)with different interatomic potentials applied at T=100 K.For COMB potential T=50 K

圖3 基于不同作用勢原子之間的距離ρ(d)和原子與其最近鄰原子之間的角度ρ(θ)的分布函數(shù)(T=4000 K;Tersoff勢T=1800 K;DL勢T=3000 K)Fig.3 Distribution functions of the distances between atoms and the angles between the lines linking one atom with its nearest neighbor ρ(d) (a)andρ(θ) (b)with different interatomic potentials applied at T=4000 K.for Tersoff potential T=1800 K and for DL potential T=3000 K

圖4顯示了在T=100 K 時基于不同原子間相互作用勢模擬所得到的原子排列。在圖4(h)中,S1和S2原子分別由黑色和綠色點表示。圖4(a)～(g)表明,具有Tersoff 勢的系統(tǒng)始終是結晶的,具有MEAM、MEAMS、SW、EDI、TM 和COMB 勢的系統(tǒng)顯然是無序的。能量的異常變化可能意味著從一種液相轉變到另一種非結晶相,如圖1 所示。應用EDI、MEAM、MEAMS、SW 和COMB 勢模擬時結果也與TM 勢類似(見圖2、3和4)。因此,從圖1、2、3和4可以得出結論,使用MEAM、MEAMS、SW、EDI、TM 和COMB勢,不能在分子動力學模擬中獲得Si金剛石結構,也就沒有晶格常數(shù)。

圖4 基于不同作用勢的原子排列(T=100 K;COMB勢T=50 K)Fig.4 Atomic configurations with different interatomic potentials applied at T=100 K.For COMB potential T=50 K(a)EDI; (b)MEAM; (c)MEAMS; (d)SW; (e)Tersoff; (f)TM; (g)COMB; (h)DL

此外,在圖2中,在DL勢的情況下,即使ρ(d)的第七個峰顯示出很小的偏差,ρ(d)和ρ(θ)也與參考的理想金剛石結構的數(shù)據(jù)一致。圖3表明,在T=3 000 K時的ρ(d)和ρ(θ)與MEAM、MEAMS、SW、EDI、TM和COMB勢相似。圖4(h)表明,每個綠色原子擁有四個黑色最近原子,并且每個黑色原子也擁有四個綠色最近原子。然而,計算結果表明纖鋅礦結構具有與金剛石結構相同的ρ(d)和ρ(θ),因此無法從ρ(d)和ρ(θ)分辨出確切的晶體結構。識別晶體結構需要檢查子系統(tǒng)的晶格,對于金剛石結構,子系統(tǒng)顯示為fcc晶格,對于纖鋅礦結構,子系統(tǒng)顯示為hcp晶格。

3.2 金剛石結構與fcc子晶格或hcp子晶格之間的關系

之前的研究表明,LJ 勢的基態(tài)對應于hcp 晶格[24]。但是,由于fcc晶格的能量比hcp晶格的能量稍高[22-23],因此在LJ勢下,系統(tǒng)還可以在某些模擬條件下顯示fcc晶格[24]。在本研究的工作中,通過更改系統(tǒng)的冷卻溫度步長n,研究n對系統(tǒng)晶體結構的影響。圖5為DL勢的分子動力學模擬中,在不同溫度步長n下,每個原子的能量和體積與溫度的關系。當n=5 K時,系統(tǒng)顯示出最低的能量和體積,以及最高的結晶溫度(1 780 K)。體積隨溫度步長的增加略有增加。

圖5 基于DL勢單原子的能量(a)和體積(b)與溫度步長n 的關系Fig.5 Dependence of energy per atom(a)and volume per atom(b)on the temperature with different temperature steps n for DL potential

圖6為在低溫時基于DL勢的分子動力學模擬中具有不同溫度步長n時的ρ(d)和ρ(θ)分布函數(shù)。圖6(a)、(c)和(e)中,ad和af分別代表金剛石結構和fcc晶格的晶格常數(shù),如箭頭所示。在圖6(c)和(e)中,對于S1(或S2)子系統(tǒng),晶胞的晶格常數(shù)為af,它對應于d2值,在該值處子系統(tǒng)在ρ(d)函數(shù)顯示第二個峰值。d1值是原子之間的最短距離。從金剛石結構的ρ(d)函數(shù)來看,金剛石結構晶胞的晶格常數(shù)ad等于af。一個金剛石晶胞中有八個原子,因此,d2和d4分別是子系統(tǒng)的d1和d2,而d1和d3分別是S1和S2原子之間的距離。d1和d2已在DL 勢中定義,但從模擬得出S1和S2原子之間的距離與DL勢中定義的距離略有不同。對于fcc子晶格,存在的角度分別為60°,90°和120°,對于hcp 子晶格,存在的角度則分別為60°,90°,109.5°和120°。

圖6 低溫時基于DL勢得到的原子之間的距離ρ(d)和原子與其最近鄰原子之間的角度ρ(θ)的分布函數(shù)(a)～(b)S1+S2 系統(tǒng);(c)～(d)S1系統(tǒng); (e)～(f)S2 系統(tǒng)Fig.6 Distribution functions of the distances between atoms and the angles between the lines linking one atom with its nearest neighborsρ(d)andρ(θ)with different temperature steps n applied at low temperatures for DL potential(a)-(b)S1+S2 system;(c)-(d)S1 subsystem; (e)-(f)S2 subsystem

當n=5 K 時,與其他冷卻速率相比,S1和S2子系統(tǒng)均出現(xiàn)第三個d峰和θ=109.5°峰共存,顯示出hcp晶格(圖6)。但當溫度步長增加時,S1和S2子系統(tǒng)中也有許多結構顯示fcc晶格。例如,n=40 K時,在加熱過程中,有一個fcc相向hcp相的轉變(圖5(a))。因此,并非所有模擬系統(tǒng)顯示hcp晶格,而是同時存在fcc和hcp晶格。

圖7顯示了DL勢在n=40 K 和T=50 K 時系統(tǒng)的原子排列。S1和S2的Si原子分別用黑色和綠色點表示,虛線框內展示的是ABAB…原子排列,而點線框內展示的是ABCABC…原子排列。Fcc點陣的原子排列是ABCABC…,而hcp點陣是ABAB…排列,其中A,B和C是密堆積的原子層,如圖7中的虛線圈所示。在n=100 K時,模擬系統(tǒng)主要是金剛石結構(圖5和6),晶化溫度(液相轉變?yōu)榫w相的溫度)為1 602 K,晶格常數(shù)為5.209?,與實驗值1 687 K和5.430?比較吻合[30]。

圖7 在n=40 K 和T=50 K 時基于DL勢得到的原子排列Fig.7 Atomic configurations for the system at n=40 K and T=50 K for DL potential

總之,應用DL 勢的系統(tǒng)在冷卻速率非常慢時得到纖鋅礦結構,而冷卻速率較快時得到金剛石結構。如果應用基態(tài)為fcc晶格的原子間勢,則會形成金剛石結構。

3.3 DL模型與實驗結果的比較

表2列出了在分子動力學模擬中應用不同原子間作用勢獲得的結晶溫度,晶格常數(shù)和彈性常數(shù)。結果表明,晶體結構應該是原子間相會作用勢有效性測試的第一個標準,只有基于DL 勢的模擬才能得到金剛石結構;對于其他勢,不能形成金剛石結構。另外,對于DL勢,彈性常數(shù)c11約比實驗數(shù)據(jù)大13%,并且c12=c44。與其他勢的結果對比,DL 勢模擬結果與實驗數(shù)據(jù)之間的一致性證明了DL勢的有效性。

目前,LJ勢已廣泛用于分子動力學模擬中,大多數(shù)研究人員認為LJ勢僅適用于稀有氣體或密堆積原子系統(tǒng)[22-23,30]。之前的研究表明,一個LJ勢能僅得到hcp或fcc晶格,而不能得到其他晶格結構,例如體心立方(bcc)晶格[24]。然而在不設置任何初始布拉格點陣情況下,通過對LJ勢的組合,可以獲得金剛石結構和石墨結構[25]、CsCl(NaCl)結構和鈣鈦礦(ABO3)結構[26],這表明LJ勢在構造原子間相互作用勢中具有重要意義。

Tersoff和EDI勢引入了三體勢來定義鍵角,并在分子動力學模擬中形成109.5°的鍵角(如圖2(b)所示)。但遺憾的是,引入三體勢并沒有獲得金剛石結構;同樣,COMB勢即使考慮了電荷因素,也無法得到金剛石結構。

在基于DL勢的分子動力學模擬中得到金剛石結構必須滿足三個條件:首先,必須創(chuàng)建兩種類型的原子,即使它們在物理上彼此相同;其次,這兩種類型的原子必須具有各自的子點陣,在這里S1和S2原子可以顯示hcp晶格或fcc晶格;第三,S1和S2原子之間的距離決定了所得的晶體結構。如果S1(S2)原子的最短距離是dA,S1和S2原子之間的最短距離是dAB,當dA/dAB之比在1.2～1.3的范圍內時,可得到CsCl結構;dA/dAB之比在1.3～1.6范圍內,可得到NaCl結構;dA/dAB之比在1.7～1.9之間得到金剛石結構;dA/dAB之比在1.9～2.1之間得到石墨結構。對于石墨結構,S1和S2原子可能顯示hcp或fcc晶格,分別可得到α-石墨或β-石墨,每個單層對應石墨烯結構。

有趣的是,基于LJ勢的分子動力學模擬中,即使基態(tài)對應于hcp晶格,當S1和S2原子之間的距離發(fā)生變化時,子系統(tǒng)的晶格也可以自適應以形成有序結構。對于CsCl結構,子晶格是dA/dAB=1.2 的簡單立方(sc)晶格;對于NaCl結構,子晶格是dA/dAB=1.4的fcc晶格,這意味著在基于LJ勢的分子動力學模擬中,系統(tǒng)可以自適應以形成能量上最低的有序結構。

4 結 論

本研究提出了適用于Si原子的分子動力學模擬的DL勢,并測試了DL 勢和其他作用勢的有效性。DL勢將LJ勢進行組合為兩個fcc晶格原子間勢和兩個晶格之間的原子間勢來描述晶體結構。在不設置任何初始布拉格點陣情況下模擬Si的晶化過程,計算模擬系統(tǒng)中原子之間的距離和原子與其最近鄰原子之間的角度的分布函數(shù)進行晶格識別。結果表明,只有應用DL勢的模擬系統(tǒng)才能獲得金剛石結構。通過應用DL勢的模擬得到的晶化溫度、晶格常數(shù)和彈性常數(shù)結果與實驗數(shù)據(jù)相吻合。