Fe团簇在Fe(110)表面上扩散和结构稳定性的分子动力学研究

王 飞， 张光华， 王三军, 3， 李金铭, 3， 王昶清， 贾 瑜

(1.郑州大学 物理工程学院 河南省量子功能材料国际联合实验室 河南 郑州 450001； 2.河南省医药学校计算中心 河南 开封 475001； 3.河南教育学院 物理系 河南 郑州 450046)

0 引言

表面原子和原子团簇扩散是表面生长和质量输运的一个重要过程.详细了解原子和原子团簇表面扩散机制对理解晶体生长、薄膜形成过程是非常重要的.很多科研工作者从实验和理论上以不同的角度研究了原子团簇在各种金属、半导体表面的结构稳定性和动力学扩散特性.

实验上通过场离子显微技术(FIM)可以研究原子团簇扩散的特性[1-4]，但是用这种技术很难直接观察原子团簇的扩散路径，尤其是原子团簇的扩散机制比较复杂的时候.因此，理论上研究原子团簇的扩散机制是必要的.原则上，研究原子团簇的扩散特性，基于密度泛函理论的第一性原理是比较精确也是首选的方法.只是由于计算机条件的限制，使用这种方法有时候受到一些限制.相对来说基于半经验势的分子动力学模拟结果不是很精确，与实验结果在定量上可能不一致，但是用这种方法来定性研究原子团簇的扩散特性还是非常有效的；而且有限温度下分子动力学模拟能够实现大尺度长时间跟踪每一个原子， 因此使用基于半经验势的分子动力学模拟研究原子团簇的扩散行为特性也是一种很好的选择[5-10].

迄今为止，相对于研究原子团簇在面心立方(FCC)金属表面扩散特性，原子团簇在体心立方(BCC)金属表面扩散的研究还比较少.最近，本文课题组采用分子动力学方法和嵌入原子势研究了金属Fe高指数表面的表面能[9]和单个空位在Fe(110)表面附近的扩散特性[10]，取得与实验比较一致的结果.文献[11-12] 也采用分子动力学方法和嵌入原子势研究了单个原子和空位在Fe(100)表面的扩散特性，也取得与实验基本一致的结果.Chen等[13]采用修正的嵌入原子势(MEAM)和分子动力学方法研究了小Fe原子团簇在Fe(110)表面扩散的动力学特性.但是，他们没有给出原子团簇的详细扩散机制.因此有必要采用分子动力学方法详细研究Fe原子团簇在Fe(110)表面的扩散机制.

本文采用嵌入原子势和分子动力学方法研究了小Fe团簇在Fe(110)表面上的结构稳定性和扩散特性.计算结果表明四原子和七原子团簇有较高的转变能和扩散势垒，说明这两种团簇比其他原子团簇较稳定，这两种团簇可能是薄膜生长过程中的临界晶核.本文还研究了三聚体团簇旋转和团簇扩散过程中二聚体剪切机制.

1 计算方法

本文中所有的计算都是采用Fe的嵌入原子势(EAM)，这一势的相关参数可以查询文献[11].尽管这一势函数是半经验的，它还是能给出Fe的一些基本特性，比如弹性常数、体模量、声子色散曲线以及单个空位和原子的扩散特性[10-12].

为了消除边界效应，计算采用20层表面超胞，每层有200个原子，平行表面方向采用周期性边界条件.垂直表面方向有2 nm的真空层.最下边两层原子固定在体位置不动，上边18层原子可以自由运动.计算过程中，采用Verlet算法积分原子的牛顿运动方程.计算体系选用NPT系综(粒子数N、压力P和温度T不变)，时间步长取2 fs.采用Andersen 等压算法[14]计算了不同温度下Fe的晶格常数以便于在有限温度下模拟时用.为了获得不同扩散路径的扩散势垒，文中采用微拉伸弹性带(nudged elastic band,NEB)算法[15].这种算法通过原子扩散的初始状态和末状态确定中间过渡态.扩散势垒由扩散过程中的鞍点过渡态和初始状态确定.为了计算原子团簇的结合能，首先把团簇的各个原子放在Fe(110)的4重对称位置，然后对体系在0 K下弛豫20 ps以便找到能量最低状态.

2 结果和讨论

2.1 原子团簇的结构稳定性

为了找到Fe原子团簇在Fe(110)表面的稳定位形，首先计算了各种位形原子团簇的结合能.原子团簇的结合能由(1)式给出[16]，

Enb=En-Es-n(E1-Es),

(1)

这里Enb，En，Es分别表示n原子团簇的结合能、衬底和n原子团簇的总能量和衬底的总能量.E1-Es代表单个原子的吸附能.本文计算的单个原子吸附能是-3.26 eV，与Chen等[13]报道的结果-3.52 eV基本符合.

图 1给出Fe(110)表面不同原子团簇的最低能位形及其结合能.从图中可以看出原子团簇的结合能随着团簇中原子个数的增加而减少.说明随着团簇原子个数增加，原子团簇的稳定性也在增强.文献[13]中也有类似的报道.

图1 Fe(110)表面Fe原子团簇的最低能位形及其结合能Fig.1 The lowest energy configurations and binding energy of Fe clusters on Fe(110) surface

一般情况下，增原子结合越紧密，其电子云交叠越多，相应的团簇能量也越低，也即是密堆团簇有较高的稳定性.从图1可以看出，几乎所有的团簇都是处于最密堆的位形.对于二聚体，计算给出最低能位形是二近邻而不是最密堆的一近邻位形，如图1(b)所示.究其原因可能是因为二近邻原子有较大的弛豫.计算也显示二近邻二聚体原子键长比一近邻的短，而且弛豫程度比一近邻大12%.对于六聚体，本文计算得到的最低能位形是六角形(图1f)，和文献[13]给出的低能方形位形不同.迄今为止也没有实验证明到底哪种位形比较稳定.然而，对于FCC(111)表面，原子团簇倾向于六角形而不是方形位形[5, 8].

为了更进一步比较原子团簇的相对稳定性，图2给出了各种原子团簇最低能位形与其他原子位形之间转变的转变能，图中箭头代表转变方向.从图中可以看出四聚体团簇和七聚体团簇有较高的转变能，分别是0.60 eV和0.70 eV(从左到右转变).另外，这两种团簇返回的转变能(从右向左)极低,只有0.08 eV和0.06 eV.尽管三聚体团簇和五聚体团簇也有较低的返回转变能(0.05 eV和0.07 eV)，但是它们的转变能比较低(0.34 eV和0.37 eV).这表明四聚体团簇和七聚体团簇是范数团簇，这两种团簇较其他团簇稳定.对于FCC金属，类似的范数团簇行为也已经被研究[1-4].

2.2 二聚体扩散

从图2(I)可以看出，单个原子沿着原子最密集方向([-111]或者[-1-11]方向)的扩散势垒是0.24 eV.这一结果与文献[13]给出的结果0.30 eV基本一致.图3给出二聚体在Fe(110)表面沿着原子最密集方向[-111]扩散的原子示意图.二聚体的其中一个原子先沿着[-111]方向迁移(如图3(a)所示)，此过程需要克服0.19 eV的势垒，随后另一原子也沿着这一方向迁移到另一个四重对称位置只需要0.16 eV的势垒(如图3(c)所示).总的来说二聚体沿着原子最密集方向([-111] 或者 [-1-11]方向)通过原子简单跳跃机制所需要的扩散势垒仅仅为0.19 eV.这一势垒要小于单个原子的扩散势垒0.24 eV，也即是说二聚体扩散要比单个原子扩散容易.这在以前的实验Pt原子在W(110)表面扩散[17]和Re原子在W(112)表面扩散[18]时也观察到类似的现象.

图2 最低能原子团簇位形的转变能Fig.2 The transition energy of the lowest energy configurations of Fe cluster

2.3 三聚体原子团簇的旋转机制

三聚体原子团簇是最小的可以拥有一维和二维的原子团簇.这里给出三聚体的旋转扩散机制，如图4所示.三聚体的最低能位形(如图4(a)所示)，其中一个原子跳到其紧邻四重对称位置(如图4(b)所示)，然后另一原子扩散到其紧邻四重对称位置(如图4(c)所示)，这样三聚体就完成了180°的旋转，这一过程需要克服0.34 eV的势垒.铝三聚体原子团簇在Al(111)表面的旋转机制在文献[5]中通过密度泛函理论计算和分子动力学模拟也进行过研究.

图3 二聚体在Fe(110)表面沿着原子最密集方向[-111]扩散的原子示意图Fig.3 Schematic of Fe dimmer diffuse on the Fe(110) surface along the most dense [-111] direction

图4 三聚体原子团簇的旋转扩散示意图Fig.4 Schematic of Fe trimer rotation diffusion mechanism on the Fe(110) surface

2.4 团簇扩散的双原子剪切模式

对于四聚体原子团簇扩散，图5给出平移和双原子剪切两种扩散模式.如图5(I)所示，四聚体原子团簇沿着表面原子最密集方向整体平移需要的势垒是0.50 eV.而通过双原子剪切模式扩散需要0.48 eV(如图5(II)所示)，比整体平移机制稍微低0.02 eV.关于双原子剪切模式，在FCC(100)表面，已有大量的工作[19-20].然而，对于BCC (110)表面，本文是第一次研究这种扩散模式.以前的工作[19]认为双原子剪切模式可以削弱原子团簇的扩散激活能.

对于更大的原子团簇(n>4)，其扩散过程是相当复杂的.希望双原子的剪切模式能够降低原子团簇的扩散势垒.图6给出较大原子团簇(从五聚体到八聚体)的双原子剪切模式和单个原子的简单跳跃机制的原子示意图.图中给出两种机制的扩散势垒，可以看出双原子剪切模式的扩散势垒比相应的单原子简单跳跃机制的扩散势垒要低.比如五聚体团簇的双原子剪切模式的扩散势垒为0.51 eV，而其相应的单原子简单跳跃机制的扩散势垒为0.64 eV.六聚体团簇的双原子剪切模式的扩散势垒为0.52 eV，而其相应的单原子简单跳跃机制的扩散势垒为0.63 eV.七聚体团簇的双原子剪切模式的扩散势垒为0.68 eV，而其相应的单原子简单跳跃机制的扩散势垒为0.74 eV.八聚体团簇的双原子剪切模式的扩散势垒为0.60 eV，而其相应的单原子简单跳跃机制的扩散势垒为0.70 eV.

图5 四聚体原子团簇平移机制(I)和双原子剪切模式(II)扩散示意图Fig.5 Schematic of Fe tetramer diffuse on Fe(110) surface((I): transition diffusion mechanism, (II): dimer shearing mechanism)

图6 团簇(从五聚体到八聚体)的双原子剪切模式和单个原子的简单跳跃机制的原子示意图Fig.6 Schematic of Fe cluster (from pentamer to octamer) dimmer shearing and sequential motion of single atom diffuse on Fe(110) surface

3 结论

采用Fe的嵌入原子势，通过分子动力学方法研究了Fe原子团簇(n<10)在Fe(110)表面结构稳定性和扩散特性.通过计算模拟可得出以下结论.

(i) 四聚体和七聚体比其他原子团簇较稳定.

(ii)单原子扩散、二聚体扩散和三原子的扩散势垒分别为0.19 eV，0.24 eV和0.34 eV.

(iii) 对于较大的原子团簇(n>3)，双原子剪切模式能够降低原子团簇的扩散势垒.

参考文献：

[1] Kellogg G L. Oscillatory behavior in the size dependence of cluster mobility on metal surfaces: Rh on Rh(100)[J]. Phys Rev Lett, 1994, 73 (13):1833-1836.

[2] Wang S C, Kürpick U, Ehrlich G. Surface diffusion of compact and other clusters: Irx on Ir(111)[J]. Phys Rev Lett, 1998, 81(22):4923-4926.

[3] Kyuno K, Ehrlich G. Diffusion and dissociation of platinum clusters on Pt(111)[J]. Surf Sci, 1999, 437(1/2) 29-37.

[4] Fu T Y, Hwang Y J, Tsong T T. Structure and diffusion of Pd clusters on the W(110) surface[J]. Appl Surf Sci, 2003, 219(1/2):143-148.

[5] Chang C M, Wei C M, Shen S P. Structual and dynamical behavior of Al trimer on Al(Ⅲ) surface[J]. Surf Sci, 2000, 465(1/2):65-75.

[6] Liu Q W, Sun Z H, Ning X J, et al. Systematical study of dimer diffusion on metal fcc(001) surfaces[J]. Surf Sci, 2004, 554(1):25-32.

[7] Wang C Q, Zhang Y S, Jia Y. A new Si tetramer structure on Si(001)[J]. Solid State Sci, 2009, 11(9):1661-1665.

[8] Flores J C, Aguilar B H, Coronado A M, et al. Double rotation mechanism in small Cu clusters concerted diffusion over Cu(111) surfaces[J]. Surf Sci, 2007, 601(4):931-935.

[9] 贾瑜，王晓春，姚乾凯，等. 金属Fe高指数面表面能的分子动力学模拟[J]. 郑州大学学报：理学版，2004, 36(2):34-38.

[10] Wang C Q, Yang Y X, Zhang Y S, et al. A single vacancy diffusion near a Fe(110) surface: a molecular dynamics study[J]. Comp Mater Sci, 2010, 50(2):291-294.

[11] Chamati H, Papanicolaou N I, Mishin Y, et al. Embedded-atom potential for Fe and its application to self-diffusion on Fe(100) [J]. Surf Sci, 2006, 600(9):1793-1803.

[12] Papanicolaou N I, Chamati H. Diffusion of a vacancy on Fe(100): a molecular-dynamics study[J]. Comp Mater Sci, 2009, 44(4):1366-1370.

[13] Chen D, Hu W Y, Yang J Y, et al. The dynamic diffusion behaviors of 2D small Fe clusters on a Fe(110) surface[J]. J Phys: Condens Matter, 2007, 19(44):446009.

[14] Andersen H C. Molecular dynamics simulations at constant pressure and/or temperature[J]. J Chem Phys, 1980, 72(4):2384-2393.

[15] Mills G, Jónsson H. Quantum and thermal effects in H2dissociative adsorption: evaluation of free energy barriers in multidimensional quantum systems[J]. Phys Rev Lett, 1994, 72(7):1124-1127.

[16] Wright A F, Daw M S, Fong C Y. Structures and energetics of Pt, Pd, and Ni adatom clusters on the Pt(001) surface: embedded-atom-method calculations[J]. Phys Rev B, 1990, 42 (15):9409-9419.

[17] Tsong T T. Experimental studies of the behaviour of single adsorbed atoms on solid surfaces[J]. Rep Prog Phys, 1988, 51(6):759-832.

[18] Stolt K, Graham W R, Ehrlich G J. Surface diffusion of individual atoms and dimers: Re on W(211)[J]. J Chem Phys, 1976, 65(8):3206-3222.

[19] Shi Z P, Zhang Z, Swan A K, et al. Dimer shearing as a novel mechanism for cluster diffusion and dissociation on metal(100) surfaces[J]. Phys Rev Lett, 1996, 76 (26):4927-4930.

[20] Trushin O S, Salo P, Alatalo M, et al. Atomic mechanisms of cluster diffusion on metal fcc(100) surfaces[J]. Surf Sci, 2001, 482 (1):365-369.