Aln(n ﹤10000)團簇熔化行為的分子動力學模擬

伊力亞爾·海米提，李春麗，方 萌，段海明

(新疆大學物理科學與技術學院，烏魯木齊830046)

1 引 言

團簇，作為幾個乃至幾千萬個原子、離子或分子的聚集體，具有獨特的物理化學性質［1］，在新材料、納米電子器件和納米催化劑等領域有重要的應用價值，已引起研究者們的廣泛關注［2-7］.對團簇的形成、結構、性質及其尺寸演變過程的研究已經成為一個新的重要領域.

團簇的熔化行為明顯不同于塊體(晶體)熔化，其熔點通常隨著尺寸(團簇所含原子數)的減小而降低，并且伴隨有明顯的預熔化溫度區間，在尺寸較小時還體現出明顯的尺寸效應(熔點隨尺寸非單調變化)及奇特的熔化行為(如固液共存態). 金屬團簇的熔化行為是人們尤為關注的性質之一. 如，Jarrold 小組研究發現熔化可以使鋁團簇的活性增強［8］;Mottet 等人發現通過摻雜可以改變銀團簇的熔點［9］;我們此前也探討了Co團簇(Co14及Co56)的表面熔化現象［10］. 此外，對金屬團簇熔化過程中的結構演化行為的研究也受到關注. 如，Schebarchov 等人在研究中發現鎳團簇在固液共存時在一定條件下團簇結構由正二十面體結構轉變為十面體結構［11］;Wen 等詳細研究了二十四面體結構的鉑團簇在熔化過程中的結構演化［12］;Xiao 等人研究發現Oh結構的Ag309團簇在存在錯位等缺陷時Oh結構向Ih結構的轉化溫度升高［13］.

三價金屬鋁晶體具有優良的導電導熱性能和優異的延展性，已經在社會生產實踐當中獲得了廣泛的應用. 對鋁團簇的研究也揭示出其一系列奇特性質，如極小尺寸時鋁團簇的鐵磁性［14］、的超穩定性及其作為未來團簇組裝材料的潛力［15］及優良的催化性能［16］等. 近年來，Jarrold小組對于中小尺寸Aln(n ＜200)團簇的熔化行為進行了大量的實驗研究，發現了一系列奇特的熔化行為:如鋁團簇的活性增強［8］，小尺寸(n＜100)時鋁團簇熔點隨尺寸變化的振蕩現象，較大尺寸(n ＞100)時鋁團簇比熱呈現雙峰的熔化現象等［17，18］. 之前我們對中小尺寸Aln(n ＜200)團簇的熔化行為進行了分子動力學模擬研究，得出了與實驗一致的結論(如小尺寸鋁團簇熔點與尺寸的反常依賴關系及中小尺寸鋁團簇熔化的比熱雙峰現象)［19，20］.

目前對于金屬團簇熔化行為的研究主要集中于中小尺寸，而對較大尺寸金屬團簇熔化行為的系統研究仍比較少見. 本文即采用分子動力學模擬方法結合退火技術，基于半經驗的Gupta 多體勢，從高對稱性(Ih、Oh)密堆積滿殼結構及退火結構出發，系統研究了包含10000 原子之內的鋁團簇的熔化行為.

2. 勢模型和計算方法

2.1 團簇初始結構

研究團簇的熔化行為，首先要確定團簇的初始結構. 一般是從某一低能結構(如基態結構)出發，研究團簇的升溫(熔化)行為. 通常初始結構可由模擬退火算法求得，但在團簇尺寸較大時，模擬退火算法找到團簇最低能量(基態)結構的效率很低. 本文在對若干Aln(n ＜10000)團簇熔化行為的分子動力學模擬中的初始結構的確定，除考慮退火算法所得結構外，還考慮了兩種常見的密堆積高對稱性結構:Ih(Icosahedron)結構和Oh(Truncated Octahedron)結構. 由于研究尺寸范圍較大，所以我們只對滿足上述兩類對稱性的滿殼層結構尺寸進行了模擬研究.

滿殼層Ih結構團簇所含原子數n 與殼層數k的關系可表示為［21］:

本文研究滿殼層Oh結構可以通過截斷正八面體的六個頂點得到. 當正八面體邊長nl和截斷距離ncut滿足關系式nl=2ncut+1 時，滿殼層Oh結構團簇所含原子數n 與Ih結構相同. 在本文研究Aln團簇尺寸范圍內(n ＜10000)，所模擬滿殼層Ih及Oh結構團簇尺寸序列為:13、55、147、309、 561、 923、 1415、 2057、 2869、 3871、5083、6525、8217，共計13 個不同尺寸.

2.2 勢函數

本文研究Al 團簇其原子間相互作用勢采用半經驗的Gupta 多體勢函數描述. Gupta 勢是基于緊束縛模型二階矩近似的嵌入原子勢，這種多體勢已被廣泛應用于研究金屬及合金金屬團簇的幾何結構［22，23］及動力學行為［24］. Gupta 勢可由原子間Born -Mayer 形式的相互作用排斥項和含多體效應的吸引項組成，形式如下

其中rij為i 、j 原子間距. r0是該原子的特征長度(文中選為Al 晶體最近鄰原子間距)，A 是衡量原子間排斥強度的量，B 是有效跳躍積分、通常只與原子類別有關. 本文研究Al 團簇其Gupta 勢參數為A =0.1221 eV、B =1.316 eV、p =8.612、q=2.516 及r0=2.86. 這些參數是通過擬合塊體Al材料的束縛能、晶格常數、體彈模量及彈性扭轉常數而得到的［25］.

2.3 分子動力學

本文采用恒溫分子動力學方法，溫度由Berendsen 熱浴標度［26］，采用速度Verlet 算法積分牛頓運動方程. 體系從初始結構出發(如上Ih及Oh滿殼層結構，及退火結構，詳見下述)，積分牛頓方程時間步長為1fs，每個溫度點進行5 ×106步MD 模擬. 模擬溫度范圍取為200K 到1200K，溫度間隔選為20K.

除考慮如上兩種高對稱性(Ih及Oh)滿殼層密堆積結構外，還采用了分子動力學模擬退火方法獲取團簇的初始結構. 由退火方法獲得團簇初始穩定結構的步驟如下:對任一團簇初始穩定結構(如可由隨機數結合最速下降法得到)直接升溫至1200K (確保團簇處于類液態)，再進行緩慢降溫的分子動力學退火過程，降溫直至200K(確保團簇處于類固態)，降溫過程溫度間隔取為20K，再對200K 時最終結構作最速下降法弛豫至穩定結構(稱為退火結構).

2.4 計算物理量

對如上每一次分子動力學模擬，在每個溫度點上計算體系的比熱. 體系比熱按下式計算:

其中Kb為Boltzmann 常數，Et為體系的總能，T為體系的溫度， ＜＞T為系綜平均(即長時間平均).

3 結果分析與討論

3.1 團簇初始結構

表1 給出13 個不同尺寸時Aln(n = 13 ～8217)團簇三種不同初始(穩定)結構的總(結合)能量. 對比分析各能量可見:在所有尺寸下Ih結構能量均為最低，Oh結構能量均高于Ih結構，考慮到鋁晶體為Oh(FCC)結構，這也說明Al 團簇在包含接近10000 個原子時其幾何結構還沒有過渡到Al 晶體，但兩類結構平均每原子結合能差值卻在減小，如從Al13的約0.0531 (eV/atom)減小到Al8217的約0.0009 (eV/atom);此外，從表1 中退火結構的能量可大致分析退火方法尋找團簇基態結構的效率問題，在團簇尺寸較小時(如309 個原子之內)，退火所得的結構就是Ih結構. 然而隨團簇尺寸增加，所得退火結構的能量均高于Ih結構，而且隨著尺寸的變大其差距也變大，表明退火結構找基態結構的效率逐漸變低，當團簇尺寸大到3871 時，退火結構的能量甚至高于Oh結構.

表1 不同尺寸時Aln (n=13 ～8217)團簇三種結構的能量(eV)Table 1 Energies (in eV)of three different structures of Aln(n=13 ～8217)clusters at different sizes

3.2 團簇的熔化行為及穩定性

圖1 (a)、(b)及(c)分別給出從三種不同初始結構(Ih、Oh及退火(Annealing)結構)出發模擬所得不同尺寸Aln(n = 309、561 及923)團簇熔化過程中的能量(E，左圖)和比熱(Cv，右圖)隨溫度(T)的變化關系圖.

對于Al309(圖1 (a))，從能量隨溫度變化圖可看出三種初始結構所對應的曲線幾無分別.這是因為Oh結構在小尺寸(309 原子)時動力學極不穩定，在很低溫度(模擬初始溫度)下已經轉變為動力學更加穩定的(類)Ih結構. 而Al309退火結構就是Ih結構(見表1 所示)，所以三條曲線都可視為是從同一個結構出發所得，三個曲線幾乎重疊. Al309的比熱隨溫度變化圖 (圖1(a)右側)出現了一個明銳的的單峰，其對應的溫度點為660K，此即為Al309團簇的熔點.

圖1 Al309 (a)、Al561 (b)及Al923 (c)團簇的能量(E)和比熱(Cv)隨溫度(T)的變化關系Fig.1 Energies (E)and the heat capacities (Cv)of Al309 (a)，Al561 (b)and Al923 (c)clusters as a function of temperature (T)

對于Al561(圖1 (b))，盡管其退火結構不是完整Ih結構(但也是類Ih結構，即相對于完整Ih結構只是部分原子有偏離)，在模擬初始溫度(200K)時，其結構也是類Ih結構，所以其曲線與Ih結構所對應曲線幾乎重合. 而Oh初始結構對應曲線在模擬初始溫度(200K)時結構已經轉變為類Ih結構，其能量曲線隨溫度的變化關系也變得與Ih結構對應曲線基本相同(為探究Al561團簇的類Oh結構向類Ih結構轉變發生時的溫度區間，我們也從更低的初始溫度(100K)出發模擬研究了Al561團簇的熔化行為，發現其結構轉變發生在160K). Al561的比熱曲線隨溫度變化圖(圖1 (b)右側)出現了一個尖銳的峰，該峰值所對應的溫度點(700K)即為Al561團簇的熔點.

由圖1 (c)可見，Al923退火結構所對應的能量曲線與Ih結構對應曲線不同(前者能量略高).由前述(表1)可知，Al923退火結構并非Ih結構，而且能量明顯高于Ih結構(平均每原子能量差也明顯大于Al561團簇的Ih結構與退火結構的能量差). Ih結構對應曲線依然保持能量最低，而Oh結構在300K 時轉變為Ih結構，這表明隨著團簇尺寸變大，Oh結構動力學穩定性變強，所以向Ih結構的轉變溫度也隨之升高. 在Al923比熱隨溫度的變化圖(圖1 (c)右側)中，300K 和720K時的兩個尖銳的峰分別對應于類Oh到類Ih結構轉變對應的溫度點及Al923團簇的熔點.

圖2 (a)及(b)分別給出從三種不同初始結構(Ih、Oh及退火(Annealing)結構)出發模擬所得不同尺寸Aln(n =2057 及8217)團簇熔化過程中的能量(E，左圖)和比熱(Cv，右圖)隨溫度(T)的變化關系圖. 從圖2 (a)可以看出，直至團簇熔化Oh結構所對應的能量隨溫度變化曲線始終都是一條獨立的曲線，這表明在完全熔化之前，Al2057始終保持了類Oh結構而沒有轉變為類Ih結構，說明Al2057團簇的(類)Oh結構在完全熔化前的全部溫度范圍內均為動力學穩定的. 在團簇尺寸n≥2057 時，(多次)模擬Aln熔化行為所得Oh結構所對應的能量隨溫度變化曲線在團簇完全熔化前都表現為獨立曲線，這表明:對Al 團簇，在包含約2000 原子以上時，其(類)Oh結構(在完全熔化發生前)是動力學穩定的.

對于其它滿殼層Aln團簇(n =2869 ～8217之間)，從三種不同初始結構 (Ih、Oh及退火(Annealing)結構)出發模擬所得團簇熔化過程中的能量和比熱隨溫度的變化關系與如上Al2057是類同的(如圖2 (b)所示Al8217情形).

圖2 Al2057 (a)及Al8217 (b)團簇的能量(E)和比熱(Cv)隨溫度(T)的變化關系Fig.2 Energies (E)and the heat capacities (Cv)of Al2057 (a)and Al8217 (b)clusters as a function of temperature (T)

圖3 五次不同模擬過程中Al1415團簇的能量(E)隨溫度(T)的變化關系Fig.3 Energies (E)of Al1415 cluster as a function of temperature (T)at five different simulations

如上所述，對于Oh結構高對稱性滿殼層Aln團簇:在所含原子數較少(n ＜1000)時，團簇動力學穩定性不高，在遠離完全熔化時就會發生結構轉變(由類Oh結構轉變為類Ih結構);而在所含原子數較多(n ＞2000)時，團簇具有高的動力學穩定性，在完全熔化發生前均不會發生結構轉變;則，對于尺寸介于1000 ～2000 之間的Aln團簇，其動力學穩定性又如何?圖3 給出由Oh初始結構出發進行5 次不同模擬所得Al1415團簇熔化過程中能量隨溫度的變化關系，作為對比，同圖中也給出了由Ih初始結構出發的相應曲線.由圖3 可見，Al1415團簇的熔化行為是非常有趣的，5 次不同模擬所得熔化行為(結構演化)不盡相同，大致可以歸為兩類情形:(1)在完全熔化前發生了結構變化，由類Oh結構轉化為類Ih結構(5 次不同模擬中有3 次屬于此類情形，見圖3 (a) ～ (c)，此3 次模擬結果的區別僅在于結構變化發生時的溫度不盡相同); (2)在完全熔化前未發生結構變化(對應于5 次不同模擬中的2 次結果，見圖3 (d)、(e)). 我們認為，這恰好體現出Oh結構Aln團簇的動力學穩定性與團簇尺寸有密切關聯:在較小尺寸(n ＜1000)時動力學不穩定，在較大尺寸(n ＞2000)時動力學穩定，而在中間尺寸(n =1000 ～2000)時，團簇動力學穩定性具有一定的隨機性，體現出一個過渡區行為.

3.3 團簇熔點隨團簇尺寸及初始結構變化

圖4 中給出了本文所研究全部Aln(n =13 ～8217)團簇從三種不同初始結構出發所得比熱曲線的峰值(團簇熔點)隨團簇尺寸n-1/3的變化關系圖，相應各熔點的具體數值列于表2 中. 由圖可見:對Al55及更大尺寸團簇，團簇熔點隨尺寸變化基本上呈現出單調變化趨勢，而在包含約200 以上原子時(Al309及更大尺寸團簇)，Al 團簇熔點與團簇尺寸的負三分之一次方呈現出近似的線性關系，該線性關系的外推值(與縱軸的交點，對應于Al 晶體情形)落在930K 附近，對應于Al 晶體的熔點(933K);由圖4 還可以看出，Al13團簇具有很高的熔點(約860K)，高于本文所研究其它所有滿殼層尺寸Aln(n =55 ～8217)團簇的熔點. 關于Al13及其附近團簇(如Al14～Al17)的高熔點以及小尺寸Aln(n =13 ～32)團簇的震蕩式的熔點尺寸變化關系，我們此前已進行了詳細的研究，得到了與實驗觀測一致的結論［19］.

圖4 Aln (n=13 ～8217)團簇的熔點隨n -1/3 的變化圖Fig.4 The melting points of Aln (n=13 ～8217)clusters as a function of n -1/3

表2 由三種不同初始結構出發時模擬所得不同尺寸Aln (n=13 ～8217)團簇的熔點(K)Table 2 The melting points (in K)of Aln (n=13 ～8217)clusters at different sizes started from three different initial structures

需要說明的是，由表2 (圖4)可以看出，對于相同尺寸Al 團簇，從不同初始結構出發時所得團簇熔點基本相同. 有些尺寸下從不同初始結構出發時的熔點雖然不相等，但是其差距沒有超出模擬所選取的最小溫度變化間隔(20K)，對此我們也對如上不同尺寸不同初始結構各Al 團簇進行了多次模擬，所得結果(相同尺寸時熔點)差別均在20K 之內，所以我們認為，對于本文所研究Al 團簇，團簇初始結構對團簇熔點沒有影響.

4 結 論

本文基于半經驗的Gupta 原子間多體相互作用勢函數，運用分子動力學方法并結合退火技術系統模擬研究了尺寸在10000 個原子之內的Al 團簇的熔化行為，分析討論了團簇能量及比熱隨溫度的變化關系. 對于各給定尺寸團簇，均考慮了三種不同的初始結構:高對稱性密堆積滿殼層的Ih、Oh結構及退火結構. 模擬結果表明:(1)在所研究團簇尺寸范圍內，二十面體(Ih)結構均顯示出高的動力學穩定性，在團簇完全熔化前無結構轉變發生;(2)Oh結構的動力學穩定性與團簇尺寸(團簇所包含原子數)密切相關:在團簇尺寸較小(1000 個原子以內)時Oh結構不穩定(熔化前會發生向類Ih結構的轉變);在中間尺寸(1000 ～2000 原子間)時Oh結構動力學穩定性發生轉變，相應于一個過渡區域;在團簇尺寸較大(2000 原子以上)時Oh結構顯示出完整的動力學穩定性; (3)對同一尺寸Al 團簇，從不同初始結構出發模擬所得團簇熔點基本相同;(4)在包含50 原子以上時，團簇熔點隨尺寸變化基本上呈現出單調變化趨勢，而在包含約200 以上原子時，Al 團簇熔點與團簇尺寸的負三分之一次方呈現出近似的線性關系.

［1］ Wang G H. Cluster Physics［M］. Shanghai:Shanghai Scientific and Technology Press，2003(in Chinese)［王廣厚. 團簇物理學［M］. 上海:上海科學技術出版社，2003］

［2］ Baletto F，Ferrando R. Structural properties of nanoclusters: Energetic， thermodynamic， and kinetic effects［J］. Rev. Mod. Phys.，2005，77:371.

［3］ Bergeron D E，Roach P J，Khanna S N. Al cluster superatoms as halogens in polyhalides and as alkaline earths in iodide salts［J］. Science，2005，307:231.

［4］ Castleman A W，Jena P. Clusters:A bridge across the disciplines of environment，materials science，and biology［J］. PNAS，2006，103:10554.

［5］ Meng K L，Zhang F S，Zhang Y P，et al. Effect of geometry structure on melting process of Na8cluster［J］.J. At. Mol. Phys.，2006，23(2):353(in Chinese)［蒙克來，張豐收，張艷萍，等.幾何結構對Na8團簇熔化過程的影響［J］. 原子與分子物理學報，2006，23(2):353］

［6］ Wang Z G，Qiu S Y，Wen Y H. Molecular dynamics simulations of the melting of nickel nanoclusters［J］.J. At. Mol. Phys.，2008，25(4):848(in Chinese)［汪志剛，邱姝穎，文玉華.納米團簇熔化過程的分子動力學模擬［J］. 原子與分子物理學報，2008，25(4):848］

［7］ Liu J T，Duan H M. Molecular dynamics simulation of structures and melting behaviours of iridium clusters with different potentials［J］. Acta Phys. Sin.，2009，58:4826(in Chinese)［劉建廷，段海明. 不同勢下銥團簇結構和熔化行為的分子動力學模擬［J］. 物理學報，2009，58:4826］

［8］ Cao B P，Starace A K，Jarrold F M，et al. Melting dramatically enhances the reactivity of Aluminum nanoclusters［J］. J. Am. Chem. Soc.，2009，131:7.

［9］ Mottet C，Rossi G，Ferrando R，et al. Single impurity effect on the melting of nanoclusters［J］. Phys. Rev.Lett.，2005，95:035501.

［10］ Lu S W，Zhang J，Duan H M. Melting behaviors of CoN(N = 13，14，38，55，56)clusters［J］. Chem.Phys.，2009，363:7.

［11］ Schebarchov D，Hendy S C. Transition from icosahedral to decahedral structure in a coexisting solid -liquid nickel cluster［J］. Phys. Rev. Lett.，2005，95:116101.

［12］ Wen Y H，Fang H，Sun S G，et al. A molecular dynamics study of shape transformation and melting of tetrahexahedral platinum nanoparticle ［J］. Chem.Rev. Lett.，2009，471:295.

［13］ Xiao X Y，Dai W C，Xia J H. The study of structural transformation for Ag309clusters with molecular dynamics method［J］. J. At. Mol. Phys.，2013，30(4):585(in Chinese)［肖緒洋，代武春，夏繼紅. Ag309團簇升溫過程中結構轉變的分子動力學模擬研究［J］. 原子與分子物理學報，2013，30(4):585］

［14］ Cox D M，Trevor D J，Whetten R L，et al. Aluminum clusters:magnetic properties［J］. J. Chem. Phys.，1986，84:4561.

［15］ Gong X G，Kumar V. Enhanced stability of magic clusters:A case study of icosahedral Al12X，X = B，Al，Ga，C，Si，Ge，Ti，As［J］. Phys. Rev. Lett.，1993，70:2078.

［16］ Shimojo F，Ohmura S，Kalia R K，et al. Molecular dynamics simulations of rapid hydrogen production from water using aluminum clusters as catalyzers［J］. Phys.Rev. Lett.，2010，104:126102.

［17］ Neal C M，Starace A K，Jarrold M F，et al. Melting of aluminum cluster cations with 31 -48 atoms:Experiment and theory［J］. J. Phys. Chem.，2007，111:17788.

［18］ Starace A K，Cao B P，Judd O H，et al. Melting of size-selected aluminum nanoclusters with 84 -128 atoms［J］. J. Chem. Phys.，2010，132:34302.

［19］ Li C L，Duan H M，Kerem M. Molecular dynamical simulation of the melting properties of Aln(n =13 -32)clusters ［J］. Acta Phys. Sin.，2013，62:193104(in Chinese)［李春麗，段海明，買力坦·開來木. Aln(n=13 -32)團簇熔化行為的分子動力學模擬研究［J］. 物理學報，2013，62:193104］［20］ Li C L，Kerem M，Duan H M. Molecular dynamic simulation of the structural and melting properties of Al196cluster［J］. J. At. Mol. Phys.，2013，30(1):91 (in Chinese)［李春麗，買力坦·開來木，段海明. Al196團簇結構及熔化行為的分子動力學模擬［J］. 原子與分子物理學報，2013，30(1):91］

［21］ Baletto F，Ferrando R，Fortunelli A，et al. Crossover among structural motifs in transition and noble -metal clusters［J］. J. Chem. Phys.，2002，116:3856.

［22］ Michaelian K，Rendón N，Garzón I L. Structure and energetics of Ni，Ag，and Au nanoclusters［J］. Phys.Rev. B，1999，60:2000.

［23］ Rodrí guez-López J L，Aguilera-Granja F，Michaelian K，et al. Structure and magnetism of cobalt clusters［J］. Phys. Rev. B，2003，67:174413.

［24］ Zhang W，Zhang F S，Zhu Z Y. Molecular dynamics study on the melting phase transition of aluminum clusters with around 55 atoms［J］. Phys. Rev. B，2006，74:033412.

［25］ Cleri F，Rosato V. Tight-binding potentials for transition metals and alloys［J］. Phys. Rev. B，1993，48:22.

［26］ Berendsen H J C，Postma J P M，van Gunsteren W F，et al. Molecular dynamics with coupling to an external bath［J］. J. Chem. Phys.，1984，81:3684.