（凝聚态物理专业论文）ga原子在gaas表面吸附行为及纳米温度计的研究.pdf

摘要 本文主要由两部分组成：a ) 单个g a 原子在g a a s ( 0 0 1 ) b2 ( 2 4 ) a s - r i e l l 表面迁移行为的研究。b ) 碳纳米管中g a 的膨胀行为的理论研究 t 在第一部分中，我们利用经验势的方法计算给出了g a 原子以及a s 原子在 g a a s ( 0 0 1 ) d 2 ( 2 4 ) a s r i c h j k 面进行迁移时的势能面。根据我们的计算， 我们找到的g a 原子在该表面上的能量较低的吸附位，这些吸附位与实验上得 到的一致。此外，我们还发现了g a 原子在该表面上进行迁移的可能的路径。 根据这些可能的迁移路径上的相对吸附能，我们发现g a 原子更喜欢沿着平行 于表面上的a s 原子对的两条路径进行迁移。对于这个结论，我们可以用c r a 原子在各迁移路径上所受的应力来解释。而a s 原予在g a a s ( 0 0 1 ) 1 3 2 ( 2 4 ) a s r i c h 表面进_ = f 亍迁移的情况则比较简单。几乎所有相对吸附能量较低的吸附 位都集中在表面a s 原子对的两侧，从而形成了平行表面a s 原子对的两条迁移 路径。另外，我们还对g a 原子在有限温度下在该表面上进行迁移的动力学过 程进行了模拟计算。计算结果与实验也是定性相符的。 在第二部分中，我们主要研究了g a 在碳纳米管中的热膨胀行为。我们首 先介绍了关于纳米温度计的概念，以及纳米温度计所使用的材料。紧接着我 们构造了c g a 之间的相互作用的势。然后将b 相的g a - 填k , 碳纳米管中。对 体系进行有限温度下的动力学模拟。所计算得到g a 在碳纳米管中的膨胀系 数，其结果与实验值接近。我们发现炭纳米管中g a 柱的长度随温度的变化基 本是线性的，这一结果有可能为纳米温度计的制备提供理论解释。 关键词：g a a s 表面经验势密度泛函理论纳米温度 a b s t r a c t t h i st h e s i sc o n t a i n st w op a 啦：a ) d i f f u s i o no fa ne x t r ag - aa t o mi nt h eg a a s ( 0 0 1 ) 陀( 2 4 ) r i e l a - a ss l l r f a c c ob ) t h ee x p a n s i o nc o e f f i c i e n t o fg r i i j u l l l e n c a p s u l a t e di nc a r b o nn a n o t u b e i np a r ta t h ep o t e - a t i a le n e r g ys u l f a c of o rt h em i g r a t i o no fa r tc x l l ag aa t o m o nt h eg a a s ( 0 0 1 ) 1 3 2 ( 2x4 ) 8 u r f l l t w a sm a p p e do u tb yp e r f o r m i n gc a l c u l a t i o n s a tt h el e v e lo f a n a l y t i c a lb o n d - o r d e rp o t e a t i a l b a s e do no u rc a l c u l a t i o n s , w eh a v e f o u n ds o m el o w e r e n e r g ys i t e sf o rt h ea d s o r p t i o no fo i le x t r ag aa t o mi nt h e s u r f a c e , w h i c ha r ei na g r e e m e n tw i t ht h ee x p e r i m e n tm o r e o v e r , m a n yp o s s i b l e p a t h w a y sf o r 趾e x t r a g aa t o md i f f u s i o ni n t h i ss u r f a c eh a v eb e e nr e v e a l e d a c c o r d i n gt ot h er e l a t i v ee l n e l g i e so f t h ep o s s i b l ep a t h w a y s ，w eh a v ep o i n t e do u t t h a tt h ei n d i v i d u a lo aa d a t o m sp r e f e r a b l yk e e pt h e i rd i f f u s i o ni nt w op a t h w a y s p a r a l l e lt ot h ea sd i m e ：l t 墨t h i sr e s u l t nb eu n d e r s t o o db yt h es t r a i nc a u s e db y t h ed i f f u s i n gg aa t o mi nt h ep a t h w a y s f u r t h e r m o r e ，w eh a v ec a l c u l a t e de x t r aa s a t o md i f f u s i o ni nc , a a s ( 0 0 1 ) s u r f l i c e ，t h er e s u l ti sm u c h s i m p l e rt i t a nt h a ti nt h e c a s eo fg aa d a t o m t w oe v i d e n tp a t h w a y sh a v eb e e nf o u n df o re x t r aa sa t o m d i f f u s i o ni nt h i ss u r f a c en e a rt h es t e p i na d d i t i o n , t h es i m u l a t e dk i n e t i cp r o c c s s e go ft h ee x t r ag aa t o md i f f u s i n gi n d i f f e r e n t p a t h w a y s a tf i n i t e t e m p e r a t u r e ss u p p o r tt h ep r e d i c t i o nf r o m o u r c a l c u l a t e dp o t e n t i a le n e r g ys u r f a c e i np a r tb ，w eh a v ei n t r o d u c e dt h en a n o m e t e rt h e r m o m e t e rt h ef i r s t l y ，a n dt h e n w ec a l c u l a t e dt h ee x p a n s i o nc o e f f i c i e n to fg a l l i u me n c a p s u l a t e di nc a r b o n n a n o t u b e w ep r o p o s e 锄a n a l y t i c a li n t e r a t o m i ep o t e n t i a lf o rm o d e l i n gt h e g a l l i u m - c a r b o ni n f r a c t i o nu s i n gas i n g l ef u n c t i o n a lf o r m t h es i m u l a t e dk i n e t i c p r o c e s s e so ft h e1 3 - c , ae n c a p s u l a t e di nc a r b o nn a n o t u b ea tf i n i t et e m p e r a t u r e sa f c c o n s i s t e n t 、析t ht h ee x p e r i m e n t t h i sr e s u l tm a y p r o v i d et h et h e o r ye x p l a n a t i o l x f o rt h en a n o m e t c rt h e r m o m e t e r p r e p a r a t i o n k e y w o r d ：g a a ss 1 1 旺a c e ，a n a l y t i c a lp o t e n t i a l ，d e n s i t y - f u n c t i o n a l t h e o r y , n a n o t h e r m o m e t e r m 第一部分g a 原子在g a a s 表面的迁移 第一章绪论 1 1g a a s 材料发展现状 化合物半导体材料砷化镓( g a a s ) 是微电子和光电子的基础材料，是化合 物半导体中最重要、用途最广泛的半导体材料，也是目前研究褥最成熟、生 产量最大的化合物半导体材料 1 】。由于砷化镓具有电子迁移率高( 是硅的 5 6 倍) 、禁带宽度大( 它为1 4 3 e v ，s i 为1 1 e v ) 且为直接带隙，容易制成 半绝缘材料( 电阻率1 0 7 1 0 9 q 册) 、本征载流子浓度低、光电特性好、以 及具有耐热、抗辐射性能好和对磁场敏感等优良特性。用砷化镓材料制作的 器件频率响应好、速度快、工作温度高，能满足集成光电子的需要 2 。 根据砷化镓材料电子迁移率高的特性，它是开发超高速计算机很理想的 器件材料。它的电子迁移率比硅高约5 倍，它的运算速度比硅器件也高许多。 在上世纪7 0 8 0 年代，人们纷纷预测并看好砷化镓材料将在超高速计算机的 发展中起着重要作用，并投入了相当的人力和财力进行研究。然而，由于一 些技术和成本问题，加之异军突起而开发出的硅材料互补型金属氧化物半导 体集成电路( c m o s ) ，其工作电压低、功耗较低、速度较快、成本低，能满足 当时的器件需要。因此，形成了对砷化镓的竞争压力，致使g a a si c 开发超 高速计算机暂时放缓了速度 3 。随着冷战结束，许多军用技术都转入民用。 由于砷化镓材料所独具的高频、高速、低噪声、低工作电压特性，它在信息 的高频高速传送和数字化的处理方面，发挥着重大的作用。用砷化镓材料所 开发出的电予器件：如金属半导体场效应晶体管( m e s f e t ) 、高迁移率晶体 管( h e m t ) 、微波单片集成电路( m m i c ) 、异质结双极晶体管( h b t ) 等在移动通 讯、光纤通讯、卫星广播通讯、情报处理以及其它一些领域发挥着硅器件不 可替代的作用，这些应用的推广大大推动了g a a s 材料的发展。 g a a s 是最重要的化合物半导体材料，它已进入快速发展时期，其晶体 生长技术和器件制作技术已较成熟，其应用领域不断扩大。g a a s 单晶材料正 向大直径、高质量、工业化大规模生产化发展。单晶直径主流产品为 7 5 m m - 1 0 0 m m ，正逐步向1 5 0 1 m 过度，直径2 0 0 m m 单晶在2 0 0 0 年初己在试验室研 究 成功。晶体质量能满足现行器件工业生产的需要。近几年，6 a a st c 增长很 快，尤其是在卫星、汽车、通信无线数据和计算机方面发展很快。估计在未 来1 0 年内，以g a a s 为主的化合物半导体材料和器件的需求量将以年均1 0 1 5 的速度增长。在该领域中，日美、欧等发达国家和地区占有绝对优势， 他们技术先进，开发和生产实力都很强，而且还在投入较多资金进一步研究 和开发。我国应随着我们的经济和国力的发展而加快该领域的发展，尽快赶 上世界水平。 1 2g a a s 材料的表面 砷化镓器件的性能除了与单晶体材料电学参数有关外还与晶片的表面 特性有很大关系。越来越多的研究表明，砷化镓的表面质量对器件的性能有 很大影响，成为制约该类器件发展的主要因素。然而在晶片加工以及随后的 器件制作过程中，几乎每步都离不开表面处理，这些处理必然会引起表面特 性的一些变化，进而影响器件的使用性能。因此。很有必要了解砷化镓晶片的 表面特性，从而有效地监控各种处理对表面的影响。 理想表面是一个假想的表面，该表面无任何杂质、缺陷，并具有沿着平 行体内晶面方向的平移对称性。最简单的理想表面在实际中很少见，因为表 面被解理以后，表面原子处于一种高度的非对称环境中，在真空一侧，由于缺 少最近邻原子，出现悬挂键。要形成稳定的表面，表面上的原子排列必须产生 畸变，以减少悬挂键的数目，于是表面原子的排列与体内原子的排列通常不 会完全相同。这种表面畸变一般表现为表面弛豫和表面重构。表面弛豫表现 为表面层原子偏离相对于理想表面原子晶格的运动，表面原胞的平移对称不 发生任何变化。最常见的表面弛豫是在表面的法线方向上产生压缩和扩张。 2 表面重构是表面原子产生偏离体晶格原子的运动后，表面层原子发生了重新 排列，表面原胞扩大。表面重构是表面化学键组合的结果。 g a a s 的表面再构的情况很复杂。在众多的再构表面中， 0 0 1 2 方向的再 构表面是最引入注目的。这是由于在大多数g a a s 材料的光电子器件中其表面 都是 0 0 1 方向的 4 。由于g a a s 是一种闪锌矿结构的半导体。故晶体的 o o t 方向是由交替的g a 原子层和a s 原子层组成，g a 原子层与a s 原子层相隔1 4 la ( 图1 1 ) 。f 0 0 1 表面是一个极性的表面，它的截面既可以是g a 也可以是a s 。 表面的最顶端的一层如果是g a 原子构成的，我们之为称作丰镓的表面 ( r i c h g a ) ：如果是a s 原子构成的，我们称为丰砷的表西( r i c h - a s ) 。当 没有再构时，截面上的g a 原子或a s 原子将形成边长为4 0a 的正方形格点。 g a a s 晶体 图1 1 0 0 1 】 s u r f a c e 4a e o o t 方向侧视图 3 o 1 3 半导体表面再构的规则 在过去的二十年中，对于半导体表面为什么会出现再构这个问题，科学 家们做了大量的理论和实验方面的工作。这些努力促使了一些关于半导体表 面再构的具有普适性的定性的规律被发现 5 - 7 。利用这些规律我们不但可 以解释大多数我们已知的表面再构的情况，还可以去预言一些未知的表面的 再构情况。关于这些规律的详细介绍在d u k e 的一本著作里给出 8 。其简单 的描述如下： p r i n c i p l ef 表面的再构总是趋向于通过重新杂化使悬挂键饱和或者把它们 转成非键态。 p r i n c i p l e2 在多数情况下，半导体性的表面要比金属性的表面稳定。 p r i n c i p l e3 表面的构型是在制备条件下能够获得的能量最低的构型。 p r i n c i p l e 4 表面趋向于自补偿 p r i n c i p l e5 表面原子的几何结构主要取决于表面原子的轨道杂化，这一杂 化是为了减少与表面键或者悬挂键相关的表面态密度。 其中规贝l j 4 常常被人称为电子计算模型( e c m ) 。这意味着一个像g a a s ( 0 0 1 ) 表面这样有极性的半导体化合物表面，它的再构必须使得悬挂键在价带上的 态是占据念，而导带上的态是空的。也就是说，表面态必须要有一个带隙来 确保表面的非金属性。值得说明的是g a a s ( 0 0 1 ) 等多数表面非常好的符合 了以上的这些规则，但是含锑或氮的半导体化合物似乎不太喜欢这些规则 9 1 2 。 下面我们来看怎样利用e h c 这个规则来判断一个给定的表面再构是否稳 定。首先要计算再构表面原胞的电子数。如果一个再构是稳定的，那么组成 这个再构的原子所能提供的电子数应与饱和表面上所有的键所需要的电子 数相吻合。 4 1 4 g a a s ( 0 0 1 ) a s - r i c h ( 2 x 4 ) 表面再构 在不同的实验条件下，g a a s ( 0 0 z ) 表面会呈现不同的再构形式。使用的 分子束外延法( m b e ) 在不同的条件下会生长出诸如r i c h a s 的c ( 4 x 4 ) ，2 4 ，c ( 2 8 ) ，1 6 ，4 6 ，r i c h - - g a 的4 2 等稳定的再构 1 3 。在这些 稳定的再构中，生长r i c h - - a s ( 2 x 4 ) 表面所需的条件是最简单，最常用的。 因此，r i c h a s ( 2 4 ) 再构也是被研究最多的g a h s ( 0 0 1 ) 表面。1 9 9 0 年， f a r r e l 和他的合作者从化学计量法的角度，对g a a s ( 0 0 1 ) a s r i c h ( 2 4 ) 表面进行了分类0 4 1 。根据他们的研究可以将a s r i c h ( 2 4 ) 表面可以分为n 相，b 相和y 相三种类型( 图1 2 ) 。 f a ) 7 ( 2 x 4 ) b ) ，( 2 x 4 ( c ，2 ( 2 x 4 ，d ) a ( 2 x 4 阑匾国匿圜匿圉 a = a t o m og 暴a t o m 图1 2 r n 1 。1 【1 i o 相是一个高温下稳定的相，它的几何结构是最上层有两个a s d i m e r ， 在a s - d i m e r 下面有一个g a - d i m e r 。8 相是一种在a n i o n r i c h 情况下稳定的构 = 饕|锨擀舷 型。它的几何结构是最上层有三个a s - d i m e r ，这是由c h a d i 最先提出的 1 5 】。 随后，n o r t h r u p 和f r o y e n 在c h a d i 的构型上稍加改进，提出了一种新的构型， 我们称为b2 1 6 。它的几何构型是最上层有两个a s - d i m e r ，在第三层有一个 a s d i m e r 。这种结构的静电能要比b 相更低，是一种比8 相更稳定的构型。 一些理论计算证明了这个结论。并且在2 0 0 0 年，l a b e l l a 等利用扫描隧道显微 镜( s t m ) 直接观察到了这个相0 7 】。y 相实际上是一个能量较高的相，它 的最顶层是个旋转了9 0 的a s d i m e r ( 图1 2 ) 。 1 5 本论文研究内容 由于b2 相的能量较低，制备比较容易。因此，b2 相是一个研究比较集 中的g a a s ( 0 0 1 ) ( 2x4 ) a s r i c h 的表面。在这些研究中，原子在表面的迁 移是一个重要的研究方向。这是由于在用m b e 生长的过程中，将不可避免的 碰到原子在表面进行迁移的情况。而对原子迁移的研究能够帮助我们更深入 的了解材料表面的基本性质。而对于这方面的理论研究并不多。 本文利用经验势结合分子动力学以及密度泛函的方法，研究了g a 原子和 a s 原子在g a a s ( 0 0 t ) p2 ( 2 x 4 ) a s - r i c h 表面的迁移行为。 6 参考文献： 1 】蒋荣华，肖顺珍，我国半导体硅材料的发展现状 j ，世界有色金属，8 ： 7 ，2 0 0 1 2 林金庭，g a a s 器件商业应用前景广阔 j ，电子材料，9 ：2 1 ，1 9 9 3 3 邵凯，一v 族化合物半导体集成电路产业发展及其应用 j ，中国集成 电路，4 4 ：3 8 ，2 0 0 2 【4 】z i a l f e r o v ，n o b e ll e c t u r c s m l ，p h y s i c s1 9 9 6 - 2 0 0 0 ，w o r l ds c i e n t i f i c p u b l i s h i n g ，s i n g a p o r e , 2 0 0 1 ，p 4 1 3 【5 】d ，j c h a d i ，t h es t r u c t u r e so f s u r f a c e si i i 【m 】，e d s s y t o n g ，m av s b h o v e ，k t a k a y a n a g ：x i e ，s p r i n g e rs e r i e si ns u r f a c es c i e n c e sv 0 1 2 4 ， s p r i n g e r , b e r l i n ( 1 9 9 1 ) ，p 5 3 2 【6 】c b d u k e , c h e m s u r f a c es c i e n c e t h ef i r s t3 0y e a r s 明，r e v 9 6 ：12 3 7 ， 1 9 9 6 7 】j p l a f e m i n a , t o t a l e n e r g yc a l c u l a t i o n so fs e m i c o n d u c t o rs u r f a c e r e c o n s t r u c t i o n s y j s u r f s c i p e p 1 6 ：1 3 3 ，1 9 9 2 【8 】c b d u k e , i np h y s i c a ls t r u c t u r e , h a n d b o o ko fs u r f a c es c i e n c e ( m 】e d i t e d b yw n u n e r t l ，v 0 1 1 ( e l s e v i e r , n o r t hh o l l a n d ) p 2 2 9 【9 】l j w l f i u n a a , p m t h i b a d o ，s c e r w i n , b r b e n n e t t , b v s h a n a b m o k , m e t a l l i ci i i - v ( 0 0 1 ) s u r f a c e s ：v i o l a t i o n so f t h ee l e c t r o nc o u n t i n gm o d e l 明， p h y s r e v l e t t 7 9 ：6 9 3 ，1 9 9 7 1 0 j e n o r t h r u p , j n e u g e b a u e r , r i m f e e n s t r a ，a r s m i r ks t r u c t u r eo f g a n ( 0 0 0 1 ) ：t h el a t e r a l l yc o n t r a c t e dc r ab i l a y e rm o d e l j ，p h y s r e v b 6 1 ：9 9 3 2 ，2 0 0 0 【l l 】a r s n f i t h , i l m f e e n s t r a , d w g r e v e ，m s s l l i i l ，m s k o w r o n s k i ，j n e u g e b a u e r , j e 。n o r t h r u p ，r e c o n s t r u c t i o n so f g a n ( 0 0 0 1 ) a n d ( 0 0 0 i ) s u r f a c e s ：o a - r i c hm e t a l l i cs t r u c t u r e s 阴，j v a c s c i t e c l m 0 1 b 1 6 ：2 2 4 2 。 1 9 9 8 【1 2 】i l m f e e n s t r a , y d o n g ，c d k ，j e n o r t h r u p ，r e c e n td e v e l o p m e n t s i n s u r f a c e s t u d i e s o f g a n a n d a l n 叨j v a c s e i t e c l m 0 1 b2 3 ：1 1 7 4 。 7 2 0 0 5 【1 3 】q 皿x u e , t ，h a s h i z u m ea n dt s a k u r a i ，s c a l 、j n 烈gt u n n e l i n g m i c r o s c o p yo fi i i - vc o m p o u n ds e m i c o n d u c t o r ( o on s u r f a c e s 阴，p r o g s u r f s e i 5 6 ：l ，1 9 9 7 1 4 】h h f a r r e l l ，c j p a l m s t r o m , r e f l e c t i o nh i 鼬e n e r g ye l e c t r o nd i f f r a c t i o n c h a r a c t e r i s t i ca b s e n c e si ng a 舾( 1 0 0 ) ( 2 x 4 卜a s ：at o o lf o rd e t e r m i n i n g t h es u r f a c es t o i c h i o m e t r y 田，j v a c s c i t e c h n 0 1 b8 ：9 0 3 ，1 9 9 0 【1 5 d j c h a d i , a t o m i cs t r u c t u r eo f o a a s ( 1 0 0 ) ( 2 1 ) a n d ( 2 4 ) r e c o r t s t r u c t e ds u r f a c e s 明， 王v a c s c i 。t e e h n o l 。a5 ：8 3 4 ，19 8 7 1 6 】j e n o r t h r u pa n ds f r o y e n , s t r u c t u r eo f g a a s ( 0 0 1 ) s u r f a c e s ：t h er o l e o f e l e c t r o s t a t i ci n t e r a c t i o n s 阴，p h y s r e v b5 0 ：2 0 1 5 ，1 9 9 4 【l7 】v p l a b e l l a , h y a n g , d w b u u o e k , p m b a d o ，p k r a t z e r ，a n dm s c h e f f l e r , a t o m i cs t r u c t u r eo f t h eg a a s ( 0 0 1 ) - ( 2 4 ) s u r f a c er e s o l v e d u s i n gs c a n n i n gt u n n e l i n gm i c r o s c o p ya n df i r s t - p r i n c i p l e st h e o r y j ， p l a y s r e v l e t t 8 3 ：2 9 8 9 ，1 9 9 9 8 第二章研究方法 在本章中我们将对本文中所采用的经验势模型、分子动力学方法以及密 度泛函理论等理论方法进行一介绍 2 1 经验势 2 1 1 经验势简介 在材料的实际生产中，远离平衡态的动态过程是不可避免的，铡如表面 的生长。研究这种动态过程，对于研究和发展材料的制作工艺有着重要的意 义。对这种动态过程的模拟般来讲有两种办法：分子动力学( 晒) 和蒙特 卡洛( m c ) 。这两种方法都是在原子尺度上对材料的动力学过程进行模拟的。 但由于计算量过于巨大，m c 和彻模拟并不合适用于一些大体系的基于第一性 原理的计算。因此，对于一些大体系来说，经验势是一种选择。经验势是指 人们通过经验的势模型，来近似地模拟原予之间的相互作用。经验势从精度 上讲，不如第一性原理那样准确，而且并不包含量子信息。但在计算量上要 远远小于第一性原理，可以快速求解大体系。如果所要研究的体系过于庞大 或者对精度的要求不高的话，经验势无疑是一个很好的选择。随着势模型的 不断改进发展，经验势的精度也越来越高，完全可以满足定性计算的需要。 科学家先后非常成功地构造了c ，s i 等半导体以及众多的金属、合金的经验 势模型 卜4 。在2 0 0 0 年，a b l e 和他的合作者提出了g a a s 的键序势 ( b o n d o r d e rp o t e n t i a l ) 模型 5 】。该势模型能够较好的研究g a a s 的各种 力学性质。这为计算g a a s 表面提供了种理论工具。 2 1 2 经验势模型 g a a s 是化合物半导体，主要通过共价键结合。由于a s 原子和g a 原子的电 9 负性不同，两者之差约为0 4 ，这导致了局域电子发生转移，从而产生有效 净电荷的相互作用。g a 单质和a s 单质则是通过混合的金属与共价键结合而成 的。因此，对于模拟g a a s 的经验势模型来讲。应能够处理金属键和共价键这 两种成键形式。在a b l e 的工作中，采用了一种短程的键序势来模拟这两种作 用。在这种势摸型中，没有考虑长程势的影响。而长程势主要来自离子与离 子之间的库仑相互作用。由于在g a h s 中，库仑势的影响远小于共价相互作用， 故长程势可以忽略。 经验势的形式如下： ：厂心)y(吩)一掣粤vae r v a ( 勺) 】( 2 - - 1 )= 厂( 勺) y ( 吩) 一= 里。；= 笠( 勺) 】 ， i j 其中e 是总能，y 8 和分别为排斥项和吸引项。 哪) = 西d oe x p - - 历( ，一) 】 嘶) = s s d 0 1 e x p 郸丽( ，_ ) 】 其中，d 。是原子对的能量，h 是原子对的平衡键长， 参数。参数b 是一个与原予对基态振动频率相关的参数。 胪后南 ( 2 2 ) ( 2 3 ) s 是一个可调节的 其中，k 是原子对基态振动的波数，u 是折合质量，c 是光速。 f ( r ) 是原子间相互用的截断函数。 ( 2 4 ) ，( ，) = r s r d j i j 1s i i l 协( r r ) ，( 2 d ) r - r d ( 2 5 ) or 屁+ d 其中，d 和r 都是可调节的参数。两原子问的距离r 若超出了r + d ，则两原 子的相互作用认为是0 。 b i 是含有角度项的键序参量。它的作用是模拟原子间的相互作用随键角 和键长的变化。公式由b r e n n e r 给出 6 ： l b o = ( 1 + 锄) 2 ( 2 - - 6 ) 锄= f ( r j ：, ) g ( o u k ) e x p 2 , u 豫( 勺一) 】( 2 _ 7 ) 七( f ，) 角度项g ( 0 ) 如下： 2 ，2c- g ( ) = 胛+ 一萨面素面丽 c 2 _ 8 ) 其中，“，c ，d 和h 都是可调的参数。 a b l e 和他的合作者分别拟合出了g a - g a ，a s a s ，g a - a s 之间的参数( 表 2 1 ) ，并且利用这些参数对各种g a 晶体，a s 晶体，g a a s 晶体进行了结合能， 晶格常数，体模，熔点等力学和热学性质的计算。计算结果与实验值、第一 性原理计算结果符合的非常好。 2 2 分子动力学方法 2 2 1 分子动力学简介 对统计力学体系进行计算机模拟时，体系在相空间中随时间的推进，需 要确定体系在各个时刻的位形。按照产生位形变化的方法，我们可以将计算 机模拟的方法分成两类： 第一类是随机模拟方法。也就是常说的蒙特卡洛方法。在此方法中，体 系位形的变化是由随机性的演化引起的。该方法可以被用到没有任何内禀动 力学模型体系的模拟上，但该方法难以处理非平衡态问题。这里我们不再详 细描述。 第二类是确定性模拟方法，即统计物理学中的所谓分子动力学方法 ( m o l e c u l a rd y n a m i c sm e t h o dm d ) 。这种方法是按该体系内部的内禀动 力学规律来计算并确定位形的转变。它首先需要建立一组分子的运动方程， 并通过直接对系统中的一个个粒子运动方程进行数值求解，得到每个时刻各 个粒子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多 体系统的静态和动态特性，从而得到系统的宏观性质。在这样的处理过程中 我们可以看出：鼢方法中不存在任何随机因素。在泐方法处理过程中方程组 的建立是通过对物理体系的微观数学描述给出。在这个微观的物理体系中。 每个粒子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论 力学上的哈密顿量或拉格朗日量来描述，也可以直接用牛顿运动方程来描 述。确定性方法是实现b o l t z m a n 的统计力学途径。这种方法可以处理与时间 有关的过程，因而可以处理非平衡态问题。 原则上，方法所适用的微观物理体系并无什么限制。这个方法适用的 体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以 是具有内部结构的体系；处理的微观客体既可以是贩子、分子，也可以是其 他微观粒子。 实际上，岫方法面临着两个基本限制：一个是有限观测时间的限制；另 一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下( 既 粒子数目趋向于无穷时) 的性质。但计算机模拟允许的体系大小要比热力学 极限小的多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们 往往引入周期性，全反射，漫反射等边界条件。当然，边界条件的引入显然 会影响体系的某些性质。 对体系的分子运动方程组采用计算机进行数值求解时，需要将运动方程 离散化为有限差分方程。数值计算的误差阶数显然取决于所采用的数值求解 方法的近似阶数。原则上，只要计算机速度足够快，内存足够多，我们可以 使计算的误差足够小。 对于仰方法，自然的系综是微正则系综，这时能量是运动常量。然而， 我们想要研究温度和压力是运动常量的系统时，系统不再是封闭的。例如当 温度为常量的系统可以认为系统放在一个热浴中。当然，在 l d 方法中我们只 是在想象中将系统放入热浴中。实际上，在模拟计算中具体所采用的做法是 对一些自由度加以约束。例如，在恒温体系的情况下，体系的平均动能是一 个不变量。这时我们可以设计一个算法，使平均动能被约束在一个给定值上。 由于这个约束，我们并不是真正的处理一个正则系综，而实际上仅仅是复制 了这个系综的位形部分。只要这个约束不破坏从一个状态到另一个状态的马 尔科夫特性，这种做法就是正确的。不过其动力学性质可能会受到这一约束 的影响。 自上世纪五十年代中期开始，如方法得到了广泛的应用。它与蒙特卡洛 方法一起已经成为计算机模拟的重要方法。应用如方法取得了许多重要的成 果，例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程 的研究。其应用己扩展到化学反应、生物学的蛋白质、重离子碰撞等广泛的 学科研究领域。 2 2 2 分子动力学模拟的基本步骤 在计算机上对分子系统的泐模拟的实际步骤可划分为四步：首先是设 定模拟所采用的模型；第二给定初始条件：第三趋于平衡的计算过程；最后 宏观物理量的计算。下面就这四个步骤分别做简单介绍。 1 模拟模型的给定 设定模型是分子动力学模拟的第一步工作。在本文中，我们采用的是上 一节中所描述的经验势模型。 当模型确定后，根据经典物理学的规律我们就可以知道在系综模拟中的 守恒量。例如对微正则系综模拟时，能量、动量、角动量均为守恒量。 2 给定初始条件 m d 模拟进入对系统微分方程组做数值求解的过程时，需要知道粒子的初 始位置和速度的数值。不同的算法要求不同的初始条件例如，g e r l e t 算法 ( 后面我们将详细介绍) 需要两组坐标来启动计算：一组为零时刻的坐标， 另一组是前进一个时间步长后的坐标，或者一组零时刻的速度值。但是，一 般来说系统的初始条件都是不可能知道的。表面上看这是一个难题。实际上。 精确选择待求系统的初始条件是没有什么意义的，因为，如果模拟足够长的 是时间，系统就会忘掉初始条件。但初始条件的合理选择，可以加快系统趋 于平衡的速度。 3 趋于平衡 按照给出的运动方程、边界条件和初始条件，就可以进行分子动力学模 拟计算。但是，这样计算出的系统不会具有所要求的系统能量，并且这个状 态本身也还不是一个平衡态。为了使系统达到平衡态，模拟中需要一个趋衡 1 4 过程。在这个过程中，我们增加或从系统中移出能量，直到系统具有所要求 的能量。然后，再对运动方程中的时间向前积分若干步，使系统持续给出确 定能量值。我们称：这时的系统已经达到了平衡。这段达到平衡所需的时间 成为弛豫时间。在如模拟中，时间步长h 的大小选择是十分重要的。它决定 了模拟所需要的时间。为了减小误差，步长h 必须取的小一些：但取的太小， 系统模拟的弛豫时间就会很长。这里需要积累一定的模拟经验，选择适当的 时间步长h 。 4 宏观物理量的计算 实际的宏观物理量往往在硒模拟的最后阶段进行。它是沿着相空间轨迹求 平均来计算得到的。例如一个宏观物理量a ，它的测量值a 。如果已知初始 位置和动量 r ( o ) 和f p 州( o ) ( 上标n 表示系综第n 个粒子的对应坐标和 动量参数) ，选择某种算法求解具有初值问题的运动学方程，便得到相空 间轨迹 ，“8 ) 和 p “ 对轨迹平均宏观物理量的表示为： 一 1 一 么= 憋亡fd f 么( ( f ) ) ，仞( f ) ) ) ( 2 - 9 ) f 呻，一，| n 、j7 如果宏观物理量动能，它的平均为 e 一： fd z e k ( pcel i m e p ( r ) ) ) = l 川( r ) ) ) t - + a o t o ( 2 1 0 ) 由于在模拟过程中计算出的动能值是在不连续的路径上的值，因此上式可表 示为在时间的各个间断点| i 上计算动能的平均值 云：击主兰业。2 州，2 ，2 一门。篇智m f 屺叫u 下面我们将讨论具体如何进行i i d 模拟。 2 , 2 3 平衡态分子动力学模拟 在经典如模拟方法的应用当中，存在着对两种系统状态的模拟。一种事 对平衡态的m d 模拟，另一种是对非平衡态的肺模拟。对平衡态系综肋模拟又 可分为如下类型：微正则系综的旧( n v e ) 模拟，正则系综的m d ( n v t ) 模拟， 等温等压系综佃( n p t ) 模拟合等焓等压系综雎dc n p h ) 模拟 7 。下面我们 仅对平衡态的m d 方法中前两类模拟做简单介绍。 一、微正则系综的m d 模拟 在进行对微正则系综的鼢模拟时，首先我们要确定所采用的势模型，这 里我们使用的是a b l e 的势模型。在这个势模型里原子问的相互作用位势是球 对称的，其哈密顿量可以写为 肚吉莩等4 - 若“饧) ( 2 - - 1 2 ) 其中，u 是势函数，是第i 个原子与第，个原子之间的距离。在这个微 正则系综中，由于这个系统的哈密顿量不显式地出现时间关联，因而系统的 能量是个守恒量。系统的体积和原子数也是不变的。此外，由于整个系统并 未运动，所以整个系统的总动量p 衡为零。这就是系统受到的四个约束。 由该系统的哈密顿量可以推导出牛顿方程形式的运动方程组 学= 土n i ， 0 升) ，。又由空间坐标可以得到原子的运动速度 形4 = ( ”+ n r i ( - - 1 ) ) 2 h ( 2 一1 7 ) 这里在n + l 步算出的速度是前一时刻的，即第n 步的速度。因而动能的计算比 势能的计算落后一步。 根据上述原理我们可以将原予数恒定、体积恒定、能量恒定的微正则系 综( n v e ) 的m d 模拟步骤设计如下： 1 给定初始空间位置 ，l ( ) ，“” ( j = l ，2 ，) 。 2 计算第n 步时原子所受的力 z ” ：一”= e 以) 。 3 利用公式( 2 1 6 ) ，计算在第n + l 步时所有原子所处的空间位置 以1 。 4 利用公式( 2 1 7 ) 计算第n 步的速度。 5 ，返回到步骤( 2 ) ，开始下一步的模拟计算。 上述算法就是所谓的v e r l e t 算法。这种形式的v e r l e t 算法动能的计算总 是比势能的计算落后一步。此外，这种算法不是自启动的。要真正求出微分 方程组( 2 1 3 ) 的解，除了需要给出初始空间位置 ，( ) 外，还要求另外给 出一组空间位置 ) 。实际上，有时候采用改进后的计算方法可能更方便： 即把n 个原子的初始位置放在网格的格点上，然后加以扰动。如果初始条件 是空间位置和速度，则采用下面的公式来计算空间的位置 ( 1 ) 0 1 = 0 妨+ 向v p + e o 乃2 2 m ( 2 一1 8 ) 然后再按上述步骤进行计算。 v e r l e t