模型二元有序合金固液界面結構的分子動力學研究*

鄭小青 楊洋 孫得彥

(華東師范大學物理系，上海 200241)

(2012年4月6日收到;2012年7月26日收到修改稿)

1 引言

固液界面結構的研究對認識凝聚態物理中的許多現象有重要的科學意義，例如固化、外延生長、浸潤、晶體生長等[1-3].當液體附著在晶體表面時，晶體的周期性結構通常會驅使其附近的液體產生有序結構.理論研究已預測了這種有序行為的存在[4-7]，其中一些預測已得到了實驗的證實[8-13].晶體引起的界面附近液體的有序結構可分為兩種類型：即垂直于界面方向的密度有序和平行于界面方向的原子二維有序.在基于X射線散射以及電子顯微技術的實驗中已經觀察到了這種垂直于界面方向的密度有序[14-17].另外，界面層液體的二維有序也在一些實驗和模擬研究中得到了證實[17].例如，在Pb和Sn的單原子層附著在Ge(111)晶體表面的簡單體系中，就能觀察到液體的二維有序結構[8,9,18].利用透射電子顯微技術，人們在Si表面吸附的液體Al-Si合金體系中也觀察到液體的二維有序化現象[19].最近人們利用高分辨透射電子顯微鏡在氧化鋁和液體鋁組成的界面體系中[20]同時觀察到液體在垂直界面方向的密度有序和平行于界面方向的二維結構有序.該實驗還發現液體的二維有序只在界面附近的三層液體中存在，明顯比垂直于界面方向的密度有序衰減得更快.

迄今為止，對界面處液體有序性的實驗研究主要集中在異質固液界面體系，即固相和液相有著完全不同的化學組成.對于固液相成分相同的界面，實驗上很難直接測量其結構特征.目前，對這類界面結構的研究主要是通過計算機模擬來實現[17,21].固液界面的模擬研究已在許多體系中開展，比如硬球體系[7]，LJ(Lennard-Jones)體系[22,23]和許多簡單金屬體系[24-26]，也包括大量對異質界面體系的模擬研究[27-30].盡管合金體系通常會更有趣，但到目前為止，對多元合金體系界面有序性的研究卻很少[31-38].

本文采用分子動力學方法研究了LJ模型二元有序合金固液界面的結構與原子擴散.選擇此體系作為研究對象，主要基于三方面的考慮.第一，對有序合金固液界面而言，界面附近通常有多個序參量同時發生變化，比如結構序、化學序(不同化學元素的規則排布)、磁有序等.這些晶體中的序參量是如何協同影響液體結構的，目前還少有報道.在當前的體系中，晶體同時具有結構和化學序的變化，有可能引起液體里結構和化學序的反常存在，這種現象在前人的研究中幾乎沒討論過.第二，該模型二元有序合金體系的固液界面是光滑界面(faceted interface).已有的對界面液體有序性的實驗或者模擬研究中，涉及的界面多是粗糙界面(rough interface)，而對合金體系光滑界面的研究卻很少.然而，對于合金光滑界面來說，界面附近液體層的二維有序結構可能更為復雜.第三，該二元有序合金代表了一大類二元有序半導體的結構.研究當前體系對認識該類半導體沿(111)方向的固化與生長有積極的意義.本文中，為獲得體系界面附近液體的結構特征，計算了垂直于界面方向的密度和結構序分布.通過分析二維徑向分布函數、二維密度和結構因子等，研究了界面層的二維結構信息.最后通過計算擴散系數研究了界面附近原子的擴散行為.我們的研究結果表明固體的結構序和化學有序對界面結構有強烈的調制作用.

2 研究方法

2.1 勢函數的選取和分子動力學模擬細節

模型二元合金由A，B兩種元素組成，它們之間的相互作用勢用所謂截斷LJ勢來描述，其中的LJ勢參數取自文獻[39］.根據該文獻，A原子與B原子之間的相互作用滿足關系：σAA=σBB，εAA=εBB且εAB=eεAA，σAB=sσAA，其中參數e=1 和s=1.2分別描述原子A與B之間的相對相互作用強度和相對距離.為了方便計算，使用一個真實的單位(氬的LJ參數)來代替約化單位.取粒子之間的相互作用勢參數εAA/kB=114.1 K和σAA=0.3405 nm.根據文獻[39］，此時體系的固液轉變溫度為72 K.有趣的是，該體系屬于congruent-melting，即在熔點液體和固體有相同的化學成分.對應于上面的勢參數，二元合金具有所謂的fcc(111)-layered結構，此結構指的是A，B兩種原子在fcc(111)方向上交替排列.該結構可以看作扭曲的閃鋅礦結構，其中(111)方向的層間距有所變化.該結構合金在二元半導體中廣泛存在，在核輻射檢測技術中所使用的melt-grown化合物半導體，例如CdTe，即屬于該類有序合金[38].在這類合金中，AB原子比較相近，這正是本文中的模型所具有的獨特性.在一定程度上，盡管簡單的LJ勢也許難以和復雜的多體勢相比，但它的形式簡單、物理明確的特點仍使它被廣泛采用.事實上，在對固液界面的研究中，簡單LJ勢往往可以獲得和實驗可比的結果[22,23].

所有的分子動力學模擬都采用LAMMPS(大尺度原子/分子的大規模并行計算)程序，時間步長取為1 fs[40].模擬體系的z軸方向垂直于固液界面，另兩個平行于固液界面的方向為x和y坐標軸，三個坐標軸方向上均采用周期性邊界條件.模擬體系的具體尺寸為76.35?A×79.34?A×338.28?A，共包含36000個原子.

為建立平衡固液共存體系，首先按照上面的尺寸產生完全有序的fcc(111)-layered結構的晶體.完整的晶體建立以后，沿著模擬體系最長方向(z方向)將大約1/2的晶體在一個非常高的溫度下熔化，而剩余的原子被固定在晶格位上.接下來，允許所有原子運動，并將整個體系的溫度降到72 K進行熱力學平衡，此過程用的是NPzAT系綜(即平行于界面方向的周期長度保持一個固定值，垂直于界面方向的周期長度是動力學變化的，為了使體系在此方向上保持零壓，該系綜是固液界面研究中常用的[27,41,42]).最后將體系在NVE系綜下進行長時間的平衡.演化平衡后體系穩定在72 K，此值即為固液共存的溫度，也即熔點.該值與文獻[39］中給出的72 K一致.在體系達到完全平衡后，在72 K體系持續50 ns的時間演化.所有的計算都是在這50 ns的平衡過程中進行的.最后建立的固液共存體系的原子坐標快照如圖1所示，圖中體系已處在固液兩相平衡態.

圖1 模型二元有序合金體系固液共存原子坐標快照 圖中紅色的實心點代表A原子，藍色的空心點代表B原子，可以看出清晰的固液界面;箭頭所指的層為界面過渡層

2.2 數密度與序參量的計算

為了描述界面的結構特征，我們計算了粒子數密度和結構序沿z方向的分布.這些分布可以刻畫對應的物理量在穿過固液界面時的變化規律.我們具體計算了兩種類型的分布：精細分布( fine-scale pro file)和經濾波法光滑的粗粒分布(coarse-grained pro files).在精細分布的計算中，首先將體系沿著z方向劃分為很多等厚度的薄層，然后計算每個薄層內的物理量的平均值，最后這些平均值隨z的變化就構成了對應物理量的精細分布.在計算中，薄層的厚度Δz取為界面法向方向層間距的1/40.為了得到各物理量的粗粒分布，采用了有限脈沖響應濾波法 ( finite-impulse-response filter，FIR)，有關該方法的技術細節見文獻[7,42].FIR方法的關鍵是選擇合適的N和ε的值，它們是控制平均加權范圍的參數.用文獻[42］的方法，對當前體系我們發現N=60和ε=145是最佳值.

數密度分布是分析固液界面處法向有序的重要物理量[32-34,42].我們分別計算了每個薄層中A原子數密度，B原子數密度與總的數密度.數密度的計算是通過50 ns的分子動力學模擬得到的.

為了分析結構變化的有序程度，還計算了結構序參量的分布.本文采用Hoyt等[25]在2001年提出的結構序參量來區分固液相原子.以前的研究表明，該序參量對bcc和fcc結構的晶體非常有效[43].根據文獻[25］，第i個原子的序參量定義為，其中，N是原子的最近鄰數目，rij是i原子與原子j之間的距離，而rideal代表該兩原子在完美晶格中的最近鄰位置.對當前關心的結構而言，取N=12.為求結構序參量的精細分布，對每個薄層里所有原子的序參量進行平均，平均歷經50 ns的分子動力學模擬.

2.3 界面二維結構的描述

除了計算各物理量沿著z方向的分布外，還研究了界面處液體在平行于界面方向的二維結構.具體地，計算了二維徑向分布函數、二維密度和二維結構因子[44].這些量的計算主要集中在界面附近的幾層原子，具體計算方法與傳統方法類似.惟一不同的是，此時原子位置僅由x和y方向的坐標決定.也就是說，此時這些界面層中的原子被投影到xy平面上.例如，二維徑向分布函數的計算與三維情況相同，除了原子i與j之間的距離被定義為最后對這一層內的所有原子求平均值，平均用的分子動力學模擬時間為50 ns.為了計算二維密度ρxy(r)，將每一層劃分為許多二維的格子(大約是200×200個)，計算每個格子里原子的平均密度.二維密度的計算通過500 ps的平均獲得.得到二維密度后，二維結構因子Fxy(k)通過下式很容易計算：

其中，ρxy(k)是每個構型的二維密度的傅里葉變換.

2.4 擴散系數的計算

為了計算擴散系數在界面附近的變化，首先將體系沿著z方向劃分為較大的區間，區間的大小大致相當于晶體中的面間距，這些區間由精細密度分布的最小值確定.對每個區間，計算該區間內每個粒子的均方位移粒子最初處在此區間中的時刻記為t=t0，每個區間的擴散系數由原子的均方位移隨時間變化的斜率得到，即：

〈〉z表示對不同時間起點(t0)的平均，這里這個平均進行了10次.為研究擴散的各向異性，本文還計算了x，y和z三個方向的擴散系數.具體計算某一方向的擴散系數時，均方位移(MSD)僅包括該方向上原子位移的貢獻.上面有關計算界面處原子擴散系數的方法在文獻[27，32—34］中常用.

3 計算結果與討論

3.1 密度和結構序的分布

圖2是精細數密度沿界面法向方向的分布.從圖中可以看出，在遠離固液界面處，三種數密度在固體區域呈規則的周期振蕩，并且A原子層和B原子層交替出現，兩種元素的密度呈現出的這種交替排列的特征是晶體結構的反映.而在遠離界面的液體區域，密度呈現平滑的分布，這正是典型合金液體的結構特征.

密度分布在界面附近的液體區域仍有明顯的波動，并且這些波動在進入液體區域大約30?A后才基本消失.這些密度波動的存在說明界面處的液體具有一定的有序性.雖然密度波動在通常的固液界面處廣泛存在，但目前的密度波動呈現出一些異質界面和單質固液界面所沒有的新特征.這些新特征體現在：第一，盡管總密度(圖2中的實線)有類似單質中的周期密度波動，但液體中的密度峰明顯不對稱，例如，圖2中d層峰的左邊明顯有“肩膀”存在，這說明界面處液體不是簡單的無序合金液體;第二，液體沒有呈現出簡單的A原子和B原子密度交替振蕩，從c層開始，總密度的波動與A，B原子密度的波動也不完全一致，例如在界面附近連續的兩層即d層和e層中，A原子的密度占主導，而c層以B原子為主;第三，微觀看來，界面處的液體層(c—g)由兩個沿z方向相互分開的亞層組成.例如，d層的總密度就是由兩個明顯分開的A原子密度峰和B原子密度峰組成.

圖2 模型二元有序合金體系固液界面密度的精細分布 實線代表總密度，點線代表A原子的密度，虛線代表B原子的密度;界面附近的幾層用 a，b，c，d，e，f和 g 標記，c 層為過渡層，對應圖1中箭頭所指的層

就異質或單質固液界面而言，液體中的密度波動主要是由固體中的周期結構引起的.而目前的系統中，固體不僅具有結構上的周期性，同時具有所謂的化學序，這正是當前界面處液體密度變化的主要原因.固體的周期結構驅使界面處液體密度的周期波動，這可由圖2中總密度的變化看出.而固體中的化學有序同樣有驅使界面處液體保留化學序的趨勢，這可以從界面附近液體A原子和B原子密度并不是完全無序看出.當然液體本身具有結構和化學無序的趨勢，這種無序的動力來自無序導致的熵增加，從而降低液體的自由能.這種固體驅使的有序和熵驅使的無序的競爭導致界面處液體的復雜的密度波動.可以看出，結構序顯然比化學序對密度的影響要強一些.因為總密度的波動和固體中的密度波動大體一致，而A原子B原子的密度波動并不完全符合固體中的密度波動.前者主要由結構序引起，而后者主要是化學序導致.這些競爭還導致一些亞層的出現，例如層c，d和e中出現的亞層對應圖2中這些層包含的兩種原子形成的兩個峰.我們還將在后面有關的結構分析中進一步闡述.

界面的有效寬度是刻畫界面特征的一個重要物理量.界面的空間分布范圍通常由粗粒密度分布的10—90寬度來衡量[33,34,42].10—90寬度定義為：塊體在固相的數密度ρs相對于液相的數密度ρl從10%變化到90%時，對應的z方向的距離.用雙曲正切函數擬合界面附近的粗粒密度分布，就可得到界面的10—90寬度.根據文獻[42］，可以計算兩種界面寬度：即固有寬度和平均寬度.固有寬度是每個時刻密度的粗粒分布對應的10—90寬度的平均值.平均寬度是將所有時刻的密度求平均，獲得精細密度分布后再求密度的粗粒分布對應的10—90寬度.

對同一個體系來說，密度固有寬度和平均寬度的差異是在原子尺度上衡量界面特征的重要指標，例如光滑界面還是粗糙界面.正如Buta等[42]討論的，這兩種寬度的差異與平衡態熔化和固化漲落導致的界面凈移動有密切關系.如果凈移動可以忽略，也就是這兩種寬度差別很小時，表明此時界面是光滑界面.相反，對粗糙界面來說，這兩種寬度存在明顯的差異，也就是說界面存在不可忽略的凈移動.可通過對比這兩種寬度來確定界面是光滑界面還是粗糙界面.通過計算，我們獲得該體系固液界面的固有寬度和平均寬度分別為12.84?A和12.91?A，說明這時界面在原子尺度上是光滑的.事實上，關于光滑界面可以清楚地從圖1中看出.

圖3 模型二元有序合金的密度和結構序的粗粒分布 實線φ?(z)和虛線ρ?(z)分別代表經過轉化、歸一的結構序和數密度分布

結構序參量也常用來研究固液界面的結構特征.圖3給出結構序參量和密度的粗粒分布.圖中的數據通過轉化、歸一使得固體對應的值為0而液體對應的值為1.通過對比可以看出，結構序和密度在界面附近彼此相交，這表明在界面附近，先是數密度的變化較快，到了一定的位置后，結構序反而變化得更快.如果我們定義結構序參量的平均值為界面中心的位置，可以看出，在靠近界面中心的固體端，結構更加有序但卻有較低的固體密度.而在靠近界面中心的液體端，雖然結構已經非常無序但卻仍有較高的液體密度.

圖4 模型二元有序合金體系的二維徑向分布函數，層的符號與圖2密度中標記的一致

3.2 界面層的二維結構分布

為了分析界面處的二維結構特征，我們首先計算了圖2中標示的b到g層的二維徑向分布函數.圖4是這些界面層對應的二維徑向分布函數gxy(r)的計算結果.從圖中可以看出，b層的gxy(r)具有典型長程有序的特點，表明這一層具有較完整的晶體結構.除了峰的高度不同外，c層和b層的二維徑向分布函數非常類似，表明c層仍然以長程有序性為主，只是有序性的強度明顯降低了.c層中，gxy(r)的形狀與一個二維有序團簇在液體環境中很相似.然而，如果c層內二維有序團簇是隨機分布的，還是不能解釋gxy(r)中呈現出的長程有序.原因可能是在c層中出現的有序團簇被b層原子提供的二維周期勢所約束，從而導致這些有序團簇不僅傾向于具有和晶體同樣的平移對稱性，同時被有效“鎖定”在“正確”的晶格位上.這種現象在Si(111)光滑界面的二維徑向分布函數中也被發現[42].由于化學序和熵的競爭，雖然d，e和g層表現出明顯的液體特征，即長程無序.但它們仍然存在著復雜的微結構，例如它們的第一近鄰峰都有明顯的“肩膀”.這種復雜結構主要是由這些層具有兩種原子形成的兩個亞層引起的，這也可以從圖2中看出.

為了深入分析界面層的二維結構，我們進一步計算了時間平均的二維密度和二維結構因子，具體結果如圖5所示.從圖5中四層原子的二維坐標快照和對應的二維密度分布中可以看出，a層和b層中的原子基本有序，d層基本顯示出無序液體的特征，而c層的結構明顯介于b和d層之間，有一定有序性的存在.

圖5 模型二元有序合金體系界面層的二維結構分布 左欄：a，b，c，d這四層中原子位置的二維快照，紅色的實心點代表A原子，藍色的空心點代表B原子;中間欄：對應層的二維密度分布，圖中線性灰色區域的點代表的密度范圍是0到0.1?A-3;右欄：對應四層的二維結構因子.k矢量在x和y方向上的范圍分別為-100×2π/Lx到100×2π/Lx和-100×2π/Ly到100×2π/Ly

二維結構因子是分析層內原子結構的重要物理量.如圖5所示，從a層到d層，晶體的對稱性的強度逐漸減小.可以看出，a層中所有k值都體現出明顯的晶體對稱性.對于b層，晶體的對稱性開始減小，但其長程有序仍然與a層的結構因子相匹配.而在c層中，較小的k值仍顯示出明顯的晶體對稱性，而較大的k值對應著液體的各向同性.而d層的結構因子顯示出基本的液體特征.從圖5中還可以看出，從固體轉變成液體就發生在界面附近的c層中.也就是說，平行于界面方向液體的結構序要比垂直于界面方向的密度有序減小得更快.

為了更直觀地理解過渡層(c層)的復雜結構，我們研究了該層原子結構的動態變化.圖6給出了該界面層的二維原子位置在幾個不同時刻的快照.圖中從上到下對應c層的結構序逐漸增強.從圖中可以看出原子形成的二維晶體團簇，這些有序團簇被周圍的液體狀的原子所包圍.從c(1)到c(3)這些團簇逐漸變大.團簇中的原子數目從幾個到幾十個不等的變化.通過對比c層在不同時刻的團簇分布圖，可以看出，結構有序是隨時間變化的，即存在團簇形成和破壞的動力學過程.這一現象與前面有關該層的結構因子和徑向分布函數相一致.

圖6 模型二元有序合金體系c層在三個不同時刻的二維原子坐標快照c(1)，c(2)和c(3)對應的時刻的結構序逐漸增強，圖中紅色的點(黑色)代表固體原子，綠色的點(灰色)代表液體原子，固液原子由結構序參量區分

3.3 擴散系數的分布

圖7是體系界面處原子的擴散系數隨界面法向距離的變化，其中圖7(a)為總的擴散系數，圖7(b)為原子在三個方向上的擴散系數.圖中標記的最開始的7層對應圖2中對應的a層到g層.從圖7(a)可以看出，從固體到液體，擴散系數單調增大，并在遠離界面的液體區達到飽和.最開始的兩層(a和b層)的擴散系數幾乎為零，表明這些層是完全有序的晶體結構，這與前面分析的結構因子的結果相符合.由于空穴和位錯俘獲在界面附近的區域中，擴散過程由液體里持續的原子運動變為界面處由空穴調制的擴散，所以原子在界面附近的擴散逐漸減慢.這種界面處擴散減慢的過程在前人研究固液界面的工作中有過類似的報道[7,45].

x，y和z三個方向的擴散系數(Dx，Dy和Dz)在遠離界面處保持一致，這與液體的對稱性相一致.在界面附近的過渡層(大約6層)中，這三個擴散分量存在明顯的各向異性：Dz＜Dx，Dz＜Dy，且Dx≈Dy.界面處擴散系數的各向異性持續了大約6層.較小的Dz值可能是由界面處液體的有序結構引起，也說明這幾層中平行于界面方向的二維擴散占主導.該結果與前人的有關計算類似[46].

圖7 模型二元有序合金體系的擴散系數 (a)總擴散系數D的分布;(b)擴散系數D的三個分量Dx(圓)，Dy(三角形)和Dz(菱形)的分布;誤差代表從100個樣本計算的平均值的二階標準偏差

4 結論

采用分子動力學方法，系統研究了模型二元有序合金體系平衡固液界面的結構和原子擴散行為.為了研究界面法向方向的結構變化，計算了數密度和結構序參量沿著界面垂直方向的分布.通過分析二維徑向分布函數、二維密度和二維結構因子研究了平行于界面方向上的二維原子層的結構.通過計算界面處原子的擴散系數，探索了界面附近原子的擴散行為.研究結果表明，當前的固液界面是光滑界面，即液體和固體存在一個明顯的分界面.界面液體層中同時存在著垂直于界面方向的密度波動和化學有序，平行于界面方向存在一定的二維結構有序.這些有序現象和單質或異質固液界面中的有序明顯不同.法向方向密度呈現A原子層和B原子層的復雜波動，這種波動形式一直延伸到液體中.對界面層的二維結構分析表明，從液體到固體的轉變就發生在界面附近靠近液體的熔融層中，液體在固液轉變層中形成了隨時間變化的二維有序團簇.在界面附近，三個方向的擴散系數分量也顯示出明顯的各向異性.

[1］Liu X Y，Boek E S，Briels W J，Bennema P 1995Nature374 342

[2］Oxtoby D W 1990Nature347 725

[3］Vlieg E 2002Surf.Sci.500 458

[4］Spaepen F 1975Acta Metal.23 729

[5］Ladd A J C，Woodcock L V 1978J.Phys.C 11 3565

[6］Mori A，Manabe R，Nishioka K 1995Phys.Rev.E 51 R3831

[7］Davidchack R L，Laird B B 1998J.Chem.Phys.108 9452

[8］ de Vries S A，Goedtkindt P，Steadman P，Vlieg E 1999Phys.Rev.B 59 13301

[9］Reedijk M F，Arsic J，de Theije F K，McBride M T，Peters K F，Vlieg E 2001Phys.Rev.B 64 033403

[10］Reedijk M F，Arsic J，Hollander F F A，de Vries S A，Vlieg E 2003Phys.Rev.Lett.90 066103

[11］Arsic J，Kaminski D，Poodt P，Vlieg E 2004Phys.Rev.B 69 245406

[12］Oh S H，Kauffmann Y，Scheu C，Kaplan W D，Ruhle M 2005Science310 661

[13］Saka H，Sasaki K，Tsukimoto S，Arai S 2005J.Mater.Res.20 1629

[14］Toney M F，Howard J N，Richer J，Borges G L，Gordon J G，Melroy O R，Wiseler D G，Yee D，Sorensen L B 1994Nature368 444

[15］Huisman W J，Peters J F，Zwanenburg M J，de Vries S A，Derry T E，Abernathy D，van der Veen J F 1997Nature390 379

[16］Cheng L，Fenter P，Nagy K L，Schlegel M L，Sturchio N C 2001Phys.Rev.Lett.87 156103

[17］Kaplan W D，Kauffmann Y 2006Annu.Rev.Mater.Res.36 1

[18］Grey F，Feidenhans’l R，Pedersen J S，Nielsen M，Johnson R L 1990Phys.Rev.B 41 9519

[19］Donnelly S E，Birtcher R C，Allen C W，Morrison I，Furuya K，Song M，Mitsuishi K，Dahmen U 2002Science296 507

[20］Kauffmann Y，Oh S H，Koch C T，Hashibon A，Scheu C，Ruhle M，Kaplan W D 2011Acta Mater.59 4378

[21］Laird B B，Haymet A D J 1992Chem.Rev.92 1819

[22］Broughton J Q，Bonissent A，Abraham F F 1981J.Chem.Phys.74 4029

[23］Huitema H E A，Vlot M J，van der Eerden J P 1999J.Chem.Phys.111 4714

[24］Jesson B J，Madden P A 2000J.Chem.Phys.113 5935

[25］Hoyt J J，Asta M，Karma A 2001Phys.Rev.Lett.86 5530

[26］Becker C A，Hoyt J J，Buta D，Asta M 2007Phys.Rev.E 75 061610

[27］Palafox-Hernandez J P，Laird B B，Asta M 2011Acta Mater.59 3137

[28］Gersermans P，Gorse D，Pontikis V 2000J.Chem.Phys.113 6382

[29］Hashibon A，Adler J，Finnis M W，Kaplan W D 2002Comp.Mater.Sci.24 443

[30］Zhang X，Rice S A 2005J.Chem.Phys.123 104703

[31］Davidchack R L，Laird B B 1999Mol.Phys.97 833

[32］Davidchack R L，Laird B B 1996Phys.Rev.E 54 R5905

[33］Sibug-Aga R，Laird B B 2002J.Chem.Phys.116 3410

[34］Sibug-Aga R，Laird B B 2002Phys.Rev.B 66 144106

[35］Becker C A，Asta M，Hoyt J J，Foiles S M 2006J.Chem.Phys.124 164708

[36］Becker C A，Olmsted D L，Asta M，Hoyt J J，Foiles S M 2009Phys.Rev.B 79 054109

[37］Ramalingam H，Asta M，van de Walle A，Hoyt J J 2002Interface Sci.10 149

[38］Henager C，Morris J R 2009Phys.Rev.B 80 245309

[39］Vlot M J，van Miltenburg J C，Oonk H A J 1997J.Chem.Phys.107 10102

[40］Plimpton S J 1995J.Comput.Phys.117 1

[41］Gao Y F，Yang Y，Sun D Y，Asta M，Hoyt J J 2010J.Cryst.Growth312 3238

[42］Buta D，Asta M，Hoyt J J 2008Phys.Rev.E 78 031605

[43］Sun D Y，Asta M，Hoyt J J 2004Phys.Rev.B 69 174103

[44］Ashcroft N W，Mermin D N 1976Solid State Physics(Toronto：Thomson Learning)

[45］Broughton J Q，Gilmer G H 1986J.Chem.Phys.84 5759

[46］Yang Y，Olmsted D L，Asta M，Laird B B 2012Acta Mater.60 4960