（凝聚态物理专业论文）高压下空位浓度对金属铝弹性性质的影响.pdf

摘要 高压下空位浓度对金属铝弹性性质的影响 凝聚态物理专业 研究生曾坤蓉指导教师何林祝文军 摘要 实际金属晶体的结构不同于理想晶体的结构，理想晶体的结构是原子完全 按空间点阵有规则的排列。实际金属晶体的结构不可避免地存在各种各样的晶 体缺陷。这些晶体缺陷对金属材料的性能会产生很大的影响。点缺陷是一类非 常重要的晶体缺陷，它们的产生、运动与相互作用，以至于聚集或消失，都影 响着金属晶体的基本性质。尤其在晶体的塑性与强度、扩散以及其它结构敏感 性的问题中扮演了重要的角色，因而长期以来对金属晶体中点缺陷的研究引起 了人们广泛的兴趣，并展开了各种理论与实验工作，其中分子动力学模拟由于 它能够从原子尺度上细致地对晶体中缺陷的性质及运动规律进行直接描述，而 得到较多应用。在已有的研究中，多为运用不同的相互作用势函数对单空位和 空位对的形成能、激化能进行模拟计算。但是目前还没有不同压力下点缺陷对 金属弹性性质定量影响的报道。 本论文运用分子动力学方法，采用嵌入原子模型( e a m ) 描述原子之间的 相互作用势，在o - 1 5 g p a 压强范围，对面心立方( f c c ) 结构金属铝中空位密 度对其弹性性质的影响进行了计算模拟。为使缺陷可以在整个空间运动，采用 了周期性边界条件。模拟均在0 k 温度下进行，对于温度的控制，采用的是 n o s 6 一h o o v e r 方法，对于应力的控制则采用了p a r r i n e l l o ，r a h m a n 等压方法。主要 工作为： 1 、运用分子动力学方法对大小为3 2 4 0 4 a 3 2 4 0 4 a 3 2 4 0 4 a ，原子数为 2 0 4 8 的系统从一1 5 g p a 到2 0 g p a 应力范围内的弹性系数进行了研究，得到了弹性 系数随压强线性变化的关系。利用m u r n a g a h a n 等温状态方程对压强和体积关系 进行了拟合，得到铝在零压下的体模量。并将p - v v 。关系与相应压强范围内的 实验数据比较，得到所采用的e a m 相互作用势的可适用的压强范围。 摘要 2 、对大小分别为2 4 3 0 3 a 2 4 3 0 3 a x2 4 3 0 3 a ，2 8 3 5 4 a 2 8 3 5 4 a 2 8 3 5 4 a ，3 2 4 0 4 a 3 2 4 0 4 a 3 2 4 0 4 a ，3 6 4 5 5 a x 3 6 4 5 5 a x 3 6 4 5 5 a 的四种晶 胞中各拿走一个点阵上的原子形成空位时各体系在零温零压下的弹性系数进行 模拟计算，得到金属铝的弹性系数随空位浓度的变化规律。 3 、对大小为3 2 4 0 4 a 3 2 4 0 4 a 3 2 4 0 4 a 的晶胞中含一个空位的系统在 0 - 1 5 g p a 应力范围的弹性系数进行了模拟，并与无空位时相同大小的晶胞的弹 性系数进行比较，得到相同空位浓度下，压强对有缺陷晶体弹性性质的影响。 关键- i z q ：缺陷：空位浓度 分子动力学；弹性系数；e a m 势 i i a b s t r a c t a b s t r a c t i nr e a l i t y , i ti sd i f f i c u l tt of i n dp e r f e c tm e t a l l i cc r y s t a l s ，w h e r ea l la t o m sa r e a r r a n g e di np e r i o d i c a ll a t t i c e m o r eo rl e s s ，v a r i o u sd e f e c t se x i s ti nr e a lm e t a l s t h e s e d e f e c t sc a ns e r i o u s l ya f f e c tt h ep r o p e r t i e so fm e t a l s a m o n gt h e m ，p o i n td e f e c ti sa n i m p o r t a n tk i n do fl a t t i c ed e f e c t s t h r o u g ht h e i rg e n e r a t i o n ，m o t i o n ，a n di n t e r a c t i o n ， t h ep o i n td e f e c t sw i l la g g r e g a t eo rd i s a p p e a ri nm e t a l s ，e s p e c i a l l yp l a ya l li m p o r t a n t r o l et oa f f e c tp l a s t i c i t y , s t r e n g t h ，d i f f u s i o na n ds t r u c t u r a ls e n s i t i v ep r o p e r t i e so f m e t a l s t h i sf i e l dh a sa t t r a c t e da t t e n t i o no fm a n ys c i e n t i s t sa n da m o u n to ft h e o r e t i c a n de x p e r i m e n t a lw o r k s h a v eb e e nd o n e i nt e r m so fm e t h o d o l o g nm o l e c u l a r d y n a m i c s ( m d ) s i m u l a t i o nh a sb e e nw i d e l yu s e dd u et ot h ea d v a n t a g et oe a s i l y d e s c r i b et h ep r o p e r t i e so fd e f e c t sa ta t o m i s t i cs c a l e b u tt h e r eh a sn or e p o r ta b o u t h o wp o i n td e f e c t sa f f e c tt h ee l a s t i cp r o p e r t i e so fm e t a l su n d e rd i f f e r e n tp r e s s u r es o f a r i nt h i st h e s i s ，t h ee f f e c to fv a c a n c yc o n c e n t r a t i o no ne l a s t i cp r o p e r t i e so f a l u m i n u mu n d e rh i g hp r e s s u r e s ( f r o m 一15g p at o2 0g p a ) i ss t u d i e db ym e a n so f m ds i m u l a t i o n t h ee m b e d d e d - a t o mm o d e l ( e a m ) i se m p l o y e dt od e s c r i b e i n t e r - a t o m i ci n t e r a c t i o ni nf a c e c e n t e r e dc u b i c ( f c c ) a l u m i n u m i no r d e rt oa v o i dt h e i n f l u e n c eo fs u r f a c e e f f e c t ，t h ep e r i o d i cb o u n d a r yc o n d i t i o n si se m p l o y e di n s i m u l a t i o n ，s ot h a tt h ed e f e c tc a nm o v ei nt h ei n f i n i t es p a c e t oc o n t r o lt h ec o n s t a n l t e m p e r a t u r ea n ds t r e s s ，n o s 6 一h o o v e rm e t h o da n dp a r r i n e l l o r a h m a nm e t h o da r e u s e d ，r e s p e c t i v e l y t h em a i nc o n t e n t so f t h ei n v e s t i g a t i o na r el i s t e da sf o l l o w s ： 1 m dm e t h o dw i t he a m p o t e n t i a lw e r eu s e dt os t u d yt h ee l a s t i cc o n s t a n t 5 u n d e r p r e s s u r e sr a n g i n g f r o m 一15 g p at o2 0 g p ai nt h ec e l lw i t hs i z eo l 3 2 4 0 4 a x 3 2 4 0 4 a 3 2 4 0 4 a t h el i n e a rr e l a t i o nb e t w e e nt h ee l a s t i cc o n s t a n t sa n c p r e s s u r ew a so b t a i n e d t h ed a t aw e r ef i t t e dt om u m a g a h a ni s o t h e r m a le q u a t i o no 。 s t a t e ( e o s ) t og e tt h eb u l km o d u l u s t h ep - v or e l a t i o ni s c o m p a r e d w i t i e x p e r i m e n t a lr e s u l ta n dw ed e t e r m i n e dt h ep r e s s u r er a n g ew h e r et h ee a mp o t e n t i a j sv a l i d 1 1 a b s t r a c t 2 f o u rc e l l sw i t hd i m e n s i o n so f2 4 3 0 3 a x 2 4 3 0 3 a x 2 4 3 0 3 a ， 2 8 3 5 4 a x 2 8 3 5 4 a x 2 8 ，3 5 4 a，3 2 4 0 4 a x 3 2 4 0 4 a x 3 2 ，4 0 4 a， 3 6 4 5 5 a 3 6 4 5 5 a 3 6 4 5 5 aa r ec h o s e n e a c hc e l lh a so n ev a c a n c yi n s i d e w h i c h c o r r e s p o n dt od i f f e r e n tv a c a n c yc o n c e n t r a t i o n t h es i m u l a t i o nc o n d i t i o ni sa to ka n d 0 g p a t h er e l a t i o nb e t w e e ne l a s t i cc o n s t a n t sa n dt h ev a c a n c yc o n c e n t r a t i o nw a s t a l c u l a t e d 3 t h ee l a s t i cc o n s t a n t so ft h ec e l l ，c o n t a i n i n go n ev a c a n c y w i t hd i m e n s i o n s o f3 2 4 0 4 a x 3 2 4 0 4 a 3 2 4 0 4 au n d e rh i g hp r e s s u r e sr a n g i n gf r o m0 g p at o15 g p a w e r ec a l c u l a t e d t h ec o m p a r i s o nw a sm a d ef o rt h ec e l lw i t ho n ev a c a n c ya n dt h e p e r f e c tc e l lw i t hs a m es i z e p r e s s u r ee f f e c to nt h ee l a s t i cp r o p e r t i e so fi m p e r f e c t c r y s t a lw a so b t a i n e d k e yw o r d s ：d e f e c t ；v a c a n c yc o n c e n t r a t i o n ；m o l e c u l a rd y n a m i c s ；e l a s t i cc o n s t a n t s e a mp o t e n t i a l i v 四川师范大学学位论文独创性及使用授权声明 本人声明：所呈交学位论文，是本人在导师鱼盐：塑塞至指导下，独立进行 研究工作所取得的成果。除文中已经注明引用的内容外：本论文不含任何其他 个人或集体己经发表或撰写过的作品或成果。对本文的研究做出重要贡献的个 人和集体，均已在文中以明确方式标明。 本人承诺：已提交的学位论文电子版与论文纸本的内容一致。如因不符而 引起的学术声誉上的损失由本人自负。 本人同意所撰写学位论文的使用授权遵照学校的管理规定： 学校作为申请学位韵条件之_ ，学位论文著作权拥有者须授权所在大学拥 霄攀位论文的部分使用权，即：1 已获学位的研究生必须按学校规定提交印刷 版和善予版学位论文，可以将学位论文的全部或部分内容编入有关数据库进行 检索：2 ) 为教学和科研目的，学校可以将公开的学位论文或解密后的学位论文 作为资料在图书馆、资料室等场所或在校园网上供校内师生阅读、浏览。 论文作者签名：吒皆匆呼二巷 泖s 辱啦黾6 & 第一章引言 1 1 概述 第一章引言 铝，钠，镁等金属是凝聚态物质的最简单形式；铝是地壳中含量最丰富的 金属元素之一，通常以铝酸盐或其它含铝化合物的形式广泛存在于矿物及岩石 之中。由于在力学、热学、电学及化学等方面的优异性能，铝及合金已经成为 现代工业，特别是航空航天领域和电子工业上的支柱性材料之一。在计算材料 科学中，经常被用作改进计算方法的试验品就是铝，它在有些范围内可以看作 是计算材料科学里的“氢原子”。因此，对金属铝性能的研究有着普遍而重要的 意义。 1 2 弹性系数计算的发展 在外力作用下任何固体( 包括晶态和非晶态) 都要发生形变乃至碎裂。当 外力撤除后，固体形变消除并恢复原状的性质称为固体的弹性性质。固体形变 时，其弹性系数是固体简谐性性质的宏观表述，对于研究固体的近远程相互作 用，物态方程和相变( 尤其是晶体稳定性) 都是很重要的。目前，晶体在常压 下弹性系数的从头计算已得到了很大的发展口。o l 。而对于如何去计算高压下的弹 性系数这一值得关注的问题却很少i l l - 1 4 】。h j fj a n s e na n d a j f r e e m a n 1 1 1 采用 全电子线形缀加平面波方法和局域密度近似方法，将总能写作晶格常数的函数， 计算了石墨的弹性系数及与等方压强的关系。t o s h i a k ii i t a k aa n dt o s h i k a z u e b i s u z a k i 1 2 1 运用平面波基和广义梯度近似( g g a ) 的第一性原理赝势方法计算 了f c c 结构固态氩压强范围为1 g p a 一8 0 g p a 的密度和弹性刚度系数。p tj o c h y m a n dk p a r l i n s k i l i 副用从头计算法由应力与应变的关系导出了溴化银的弹性常数， 并得到了其随压强的变化关系。o g f i l s e r e na n dr e c o h e n h 1 应用线形缀加平 面波方法和混合基赝势方法计算了压强高达4 0 0 g p a 、温度达1 0 0 0 0 k 的体心立方 结构的钽的热弹性性质。gvs i n k o 和nas i m r n o v 1 5 1 运用从头计算法，根据应 变与系统总能的关系导出了铝的体积与弹性常数的关系。以上这些研究都采用 理想晶体结构，真实材料中的缺陷是不可避免的，对于缺陷在不同压力下对弹 第一章引言 性系数的影响还少见报道。 1 3 点缺陷 实际的晶体中总是存在着各种各样的缺陷，偏离了理想晶格的情况。最近 几十年来，固体科学技术的发展，逐步地、越来越深入地揭示出，在晶体内部 存在着各种各样的缺陷，它们对于晶体的各种性质产生着十分重要的影响。 1 3 1 点缺陷的分类、产生及对晶体一陛质的影响 原子尺度的晶格缺陷称为点缺陷，其特征是各个方向的尺寸很小，也称为 零维缺陷。晶体中的点缺陷包括空位、间隙原子、杂质原子及其组合成的复杂 缺陷【1 6 _ 1 7 】 1 ) 空位 在晶体中，位于点阵格点的原子在其平衡位置不停地做热振动，在任一瞬 时，总有些原子的能量大到足以克服周围原子对它的束缚，从而脱离原来的平 衡位置迁移到别处，在原来的位置上出现了空格点。从晶体中正常格点位置上 取走一个原子形成空位，如果此原子移动到晶体表面或晶体内部的某界面， 所产生的空位称为肖脱基缺陷。在晶体内部的原子也可以由于热涨落由格点跳 进间隙位置，从而产生个空位和个间隙原子。我们将这样的一对缺陷称为 弗伦克尔( f r e n k e l ) 缺陷，产生的空位称为弗伦克尔空位。平衡状态下空位的浓 度是温度的函数，温度升高，则原子脱离原来位置的几率增大，空位浓度提高， 反之，温度降低，空位浓度减少。 2 ) 间隙原子 如果在点阵的间隙位置挤进一个同类原子，就形成一个间隙原子，间隙原 子通常有三种构形：1 ) 填充在晶体的间隙位中，如八面体间隙和四面体间隙两 种；2 ) 两个原子共占一个晶格位，即哑铃构形( d u m b b e l lc o n f i g u r a t i o n ) ：3 ) 挤 列子构形( c r o w d i o nc o n f i g u r a t i o n ) ，即在最密排原子线上由n + 1 个原子共占n 个晶格位置。 如果在晶体中挤进一个异类原子，当异类原子的体积与基体原子相差不多 时，则异类原子以置换形式存在，若体积相对较小时，则有可能以间隙原子的 2 第一章引言 形式存在。 然而，在某些条件( 如在淬火、范性形变或辐照之后) 下，会存在大量过 饱和的点缺陷，这些点缺陷能够相互作用，当一个空位和一个间隙原子结合时， 间隙原子重新占据了正常的格点位置，这两个缺陷就会消失。如果两个空位结 合成双空位，就形成了最简单的缺陷集团( d e f e c t c l u s t e r s ) ，这个缺陷集团继续吸 收其它缺陷，进而形成小的空穴( v o i d ) 。 利用高温淬火，高能粒子辐照，合金化和塑性形变等方法都可以产生缺陷。 在面心立方( f c c ) 结构的金属中，由于形成空位的能量只是形成间隙原子的三 分之一，空位是比间隙子更容易形成的缺陷1 8 l ，本文研究的对象就是空位缺陷。 点缺陷的存在，其一系列的物理性能都会受到影响1 1 9 1 ，比较引人注意的是点缺 陷对于密度及电阻的影响： ( a ) 密度和线度如果在点阵中取走一个原子，放到表面上去，点阵就形 成一个空位。如空位周围原子都不移动，则应使晶体体积v 增加一个原子体积， 而点阵参数a 保持不变。但是实际上空位周围原子会产生位移，因此v 和a 都 有变化。理论计算的结果表明，间隙原子引起的体膨胀约在1 也原子体积，而 空位的体膨胀约为o 5 原子体积。 ( b ) 电阻点缺陷对于传导电子产生附加的散射，引起电阻的加大。关于 点缺陷的电阻的理论计算，一般是套用夫里德耳的合金理论，将点缺陷看为零 价或一价的杂质的原子。但是间隙原子所引起的畸变较大，效应不易正确估计， 结果的差异也较大见表1 1 。 表1 1 点缺陷产生的附加电阻( 微欧姆一厘米9 6 ) 1 2 0 l 第一章引言 在一般情况下，点缺陷对金属力学性能的影响较小。它只是通过和位错交 互作用，阻碍位错运动而使晶体强化。但在高能粒子辐照的情形下，由于形成 大量的点缺陷和挤列子，会引起晶体显著硬化和脆化，这种现象称为辐照硬化。 1 3 、2 点缺陷模拟的进展 空位是固体中最简单和常见的一种缺陷，但是固体的许多物理和机械性质 与空位的存在关系极大。详细地研究空位的形成和性质对理解材料的微观和宏 观性能都有重要作用。因此，高压条件下，对于铝的弹性系数的研究有广泛而 重要的意义，而真实材料都不可避免的存在各种各样的缺陷，于是对含有虽然 简单但对材料性质又有重大影响的空位的晶体在高压下的弹性性质又有着深远 的意义。 自8 0 年前，f r e n k e l l 2 t 首先对热平衡中的晶格缺陷介绍后，晶格缺陷的研究， 不管是实验的还是理论的就没有停止过。实验方面，与六七十年代研究晶体缺 陷一样，用辐照损伤2 2 和放射实验 2 3 】方法研究了许多金属系统。尽管不是很直 接，但仍然可以比照晶体中的情况，运用相同的工具指出相同缺陷存在的一切 可能性。正电子湮灭方法【2 4 - 2 5 1 在这方面也很有用。 早期，利用分子动力学研究空位形成能，空位迁移能。其空位能及迁移能 的准确性决定于相互作用势的准确程度。早期应用分子动力学对缺陷的研究中 多采用对势方法。b e n e n e t t 2 6 1 等应用分子动力学对l e n n a r d - j o n e s ( l - j ) 玻璃体 中的空位和空位团簇进行了研究，发现在l j 玻璃体中缺陷的迁移温度很低，这 些空位和空位团簇在很低的温度下就消失了。e m a n u e l ep a c i 和g i o v a n n i c i c c o t t i 【2 7j 采用分子动力学方法研究了f c c 结构的氩在温度为6 0 k 和8 0 k 时的空位 跳跃速率常数，其原子间的相互作用势采用的是l e t m a r d - j o n e s 两体势进行描述。 l n d e reb a t r a 和f a r i dea b r a h a m i 28 】用分子动力学和包含两体和三体形式的 s t i l l i n g e r - w e b e r 模型势模拟计算了硅晶体中的自间隙子的原子弛豫和弛豫能。 t o m o n o r ik i t a s h i m a 等 2 9 1 也采用了包含两体和三体形式的s t i l l i n g e r - w e b e r 模型势 研究分析了砷化镓中点缺陷的扩散，得到即使在熔点附近，固体中空位的扩散 也很微弱，自间隙砷原子的扩散系数比自间隙镓原子的值大。随着计算科学的 发展，有了更接近真实的多体势函数，其中包括嵌入原子方法( e a m ) 和 4 第一章引言 f i n n i s s i n c l a i r 类型的多体势。d a w 和b a s k e s 提出嵌入原子方法( e a m ) 之初口， 空位形成能就作为拟合势函数参数的重要物理量，并模拟了在镍晶体中形成4 个 空位所引起的层裂。后来，在此基础上他们运用分子动力学和嵌入原子方法进 一步计算得到了许多基态性质，如晶格常数、弹性系数、升华能和空位形成能： 并研究了金属镍和钯中的杂质、表面及其它缺陷【3 l 】。随后运用多体势对金属性 质进行研究中，空位的形成能及迁移能对多体势函数的参数值拟合非常的重要， 于是对其计算也不可少。b a s k e s 3 2 1 采用包含了角度力、势的修正嵌入原子势 ( m w a m ) 对m o s i 系统进行了研究，并用相稳定性，晶格常数，弹性系数和 点缺陷( 空位和间隙子) 能、面缺陷能对势函数进行了预测，结果与实验符合 的很好。 在运用第一性原理研究缺陷方面，首先取得突破的是c a r 和p a r r i n e l l o l 3 ， 他们应用总能梯度( 与k o h n s h a m 轨道相关) ，结合平面波基和规范存储离子的 赝势引入了迭代方法。随着计算机模拟技术的不断提高，利用第一性原理研究 铝的空位对其各方面性质的影响也越来越多：g i l l a n f 3 4 j 计算了简单金属的空位 形成能，著给出了铝的结果，他的计算是基于赝势方法，交换关联的局域密度 近似( l d a ) 以及超原胞几何学。rwj a n s e n 3 5 】等用赝势方法中交换关联及局 域密度近似方法计算了铝和铌的堆积性能和铝的空位形成能和高对称性正八面 体的铝的自间隙形成能。m j m e h la n db m k l e i n l 3 6 1 采用全电子线形缀加平面波 ( l a p w ) 方法来完成在超原胞近似中铝的空位形成能的自洽总能电子结构的 计算；计算了4 ，8 ，1 6 ，2 7 个点阵的超原胞未弛豫和弛豫后分别的空位形成能。 nc h e t t y l 3 7 等用第一性原理的赝势方法和大的超原胞计算了镁和铝的空位形成 能，并讨论了镁和铝含有杂质对其性能的改变。k a r i nc a r l i n g 和g 6 r a n w a h n s t r 6 m 等 3 8 - 3 9 1 采用局域密度近似( l d a ) 和广义梯度近似( g g a ) 方法精 确的描述了铝中空位，解决了因周围电子边界的关联效应产生的不足，讨论了 双空位的力学不稳定性。 在这些研究中，多是运用没有可调参数的第一性原理方法，选用不同的赝 势，应用超原胞几何学选择不同数目的晶格点阵数目，( 如朱梓忠【4 0 】就详细讨论 了空位形成能的理论计算值与所用超原胞大小，k 点数等参数的关系。) 以及不 同交换关联方法( l d a 或g g a ) 对铝中空位的形成能及相关力学性质进行了研 第一章引言 究。 从以上可见，由于真实材料不可避免的存在缺陷，人们对于缺陷研究的热 情一直以来都很高，用不同的方法对缺陷及缺陷对晶体的各种影响进行了研究， 但是对于缺陷对晶体弹性性质的研究还少见报道，特翱是对于高压下，缺陷浓 度对金属晶体的弹性性质的定量研究还不曾多见。 1 4 本文将进行的工作 本文将运用分子动力学方法对不同压力下不同空位浓度对金属铝的弹性性 质的影响进行计算机模拟。铝的相互作用势用e a m 方法描述，势函数形式和各 参量的值均采用y m i s h i n a 1 - 4 2 1 给出的结果，分别对无缺陷的铝在高压下的弹性 系数，四种含不同空位浓度的系统在零压下的弹性系数及高压下缺陷系统的弹 性系数进行模拟，得到高压下空位对弹性性质影响的定量关系。 6 第二章模型方法 2 1 分子动力学方法 第二章模型方法 分子动力学方法以其能够跟踪粒子轨迹、模拟结果准确的特性而倍受研究 者的关注，其基本原理是假定原子的运动是由牛顿运动方程决定的。将由n 个 原子组成的体系抽象成相互作用的质点，然后给出这n 个质点间的相互作用， 运用经典力学方程求解每个粒子的运动轨迹，在此基础上研究体系的结构及其 它相关性质，因此该方法适应于研究材料的结构和性质随时间或温度变化的过 程。分子动力学模拟给出系统中任一时刻的所有粒子的坐标和速度的微观信息， 由此可按统计力学平均得到系统的热力学参量。a i d e r 和w a i n w r i g h t l 4 “4 吁1 9 5 7 年使用硬球模型实现了第一个分子动力学模拟计算。1 9 6 0 年g i b s o n 等m 】使用分 子动力学探讨材料辐射损伤的动力学规律，给出了原子轨迹，这一工作使得过 去对材料性能的定性估计迈向微观结构和过程的定量研究。1 9 6 4 年r a t u n a n 【4 6 j 用m d 方法模拟液体，同时加进了周期性边界条件，结果惊喜的发现可以用很 少的粒子( 8 6 4 个) 来反映真实系统的热力学性质。自此，凝聚态物质的分子动 力学模拟成为可能，许多研究者开始了运用分子动力学方法的研究。分子动力 学模拟的主要用途有：液体模拟，晶体中的缺陷模拟，材料的断裂模拟，表面 物理，摩擦学，分子生物学和制药，材料的电学特性和动力学性质。下面对分 子动力学方法进行概述。 2 1 1 分子动力学的基本运动方程 分子动力学是一个确定性的数值模拟方法，从理论上说，只要给定各个原 子的初始位置、速度和边界条件，之后系统的时间演化就完全确定了。它是通 过一系列牛顿方程的解来很好的定义的系统的经典统计力学。具有坐标再，速 度耳和质量肌，i = 1 ，n 的粒子数为n 的经典体系，其中的粒子是通过势 ( 亏，两) 来相互作用的。大多数情况下可写作下列形式： 第二章模型方法 = “j ( 2 1 ) 其中，勺= l 弓l = 睇一弓l ，还可以有更为普搀的形式，可能会花更多的计算时 间。牛顿方程为： m z i2 篇z l 。嘞d 庐i ，f = l ， ( 2 2 ) 分子动力学的基本任务就是求解该方程的数值解，粒子的运动轨迹通常由 差分方法获得。其基本思想是把模拟时间按一个固定时间步长a f 分割。假定在 f t + 西时间段内粒子受到的力保持不变，由f 时刻粒子的速度和坐标，以及所 受的力，可以得到t + 西时刻粒子的位置和速度，新的受力状态又可由此决定， 重复循环，直到系统经历统计所必须的所有状态。感兴趣的统计平均是从i p ，和 巧一，i = 1 ，n 计算而来的。巧f f j 和i r f ，是系统在相空间整个轨迹的时间 平均。由于某些原因，n 被限制在几千的数值，所以对于小系统来说，表面效应 显然是很重噩。 2 1 2 有限差分法 预测一校正法 解像方程( 2 1 ) 和) ( 2 2 ) 这样普遍的微分方程一种标准的方法就是有限 差分法。基本思路如下，给定分子在时刻t 的位置，速度和其它动力学信息， 尝试得到在下个时刻t + 6 t 精确度足够高的位置，速度等。方程是一步一步来 求解的，时间间隔所的选取多少依赖求解的方法，不过新肯定比一个分子移 动它自身线度所需的时间小。 如果经典轨迹是连续的，那么时刻t + f i t 的位置，速度等可以由时刻t 时的相应物理量的泰勒展开得到： f 9 0 + 6 t ) = i o ) + 占f 伊o ) + d f2 丘o ) + 占f3 占( f ) + ( 2 3 ) 上 o f 9o + 占，) = 旷o ) + 占，厅0 ) + 占f 2 占0 ) + - ( 2 4 ) 第二章模型方法 厅，o + 况) = 厅o ) + 所方( f ) 十 ( 2 5 ) 占9 0 + c y t ) = 占o ) ( 2 6 ) 上标表示这些是“预测值”，马上要对它们进行“校正”。但是上式并不能 随时间的进行雨产生校正轨迹，因为并没有引入运动方程。可以从薪的位置i 来 计算时刻f + 艿r 的力和校正加速度酽o + 舀) ，可将这些值与上式的预测加速度 比较，估算出预测步骤误差的大小： 五o + 拉) = 石。( f + 国) 一a ，( t + 6 t ) ( 2 7 ) 将这个误差连同预测步骤的结果反馈到校正步骤，就有： i ( f + 国) = f 9 0 + 詹) + c 。石o + 舀) ( 2 8 ) f 。( f + 舀) = i 9 0 + 痨) + c 厅o 十痨) ( 2 9 ) 五。0 + 拉) = 厅9 ( f + 痨) + c ：厅( f 十揖) ( 2 1 0 ) 占。妇+ 函) = 5 ，o + & ) + c ，a ( f + ) ( 2 1 1 ) 于是f ( f + 西) 等物理量更好的接近了真实的位置，速度等。g e a r 4 7 。4 8 1 给出了 上述方程中系数的最佳值，截去了泰勒展式中三阶以上的导数值，这种情况的 最佳系数值为c 。= 1 6 ，c 。= 5 6 ，c ：= 1 和c ，= 1 3 。校正步骤是一个不断重复的步 骤：由位置f 。并与当前的a 。比较计算出更新的“校正”加速度，如此，通过像 方程( 2 - 8 ) 一( 2 - n ) 就可以进步使位置，速度等更精确。在很多应用中，这 样的重复是得到精确解的关键。预测值为方程的解提供了初始猜想，但它在理 论中却不是很好的解，因为连续的校正重复会很快速的收敛到正确结果。在分 子动力学中，由于每步都暗含校正步骤，大量的校正重复是最费时间的。通常 的，只有个( 偶尔两个) 校正步骤可解，所以精确的预测阶段( 就如方程( 2 8 ) 一( 2 1 1 ) 中的泰勒展开式) 是很重要的。基于预测校正法的逐步分子动力学模 第二章模型方法 拟的基本方案如下： ( a ) 用当前的位置、速度、加速度等值来预测在时刻，+ 万r 的对应各值； ( b ) 从新的位置来计算力，于是也可计算出加速度丘。= z m 。； ( c ) 用新的加速度来校正预测过的位置、速度、加速度等； ( d ) 计算感兴趣的变量，如能量、势能、量级参数，在返回步骤( a ) 进 行下一步预测校正前，为时间平均的积累做准备。 v e r l e t 算法 应用最广泛的运动方程的积分最早是由v e r l e t 【4 9 1 提出并归功于s t s r m e r d 8 1 的 v e r l e t 算法这是一种直接求解二阶微分运动方程的方法。此方法是基于位置 f ( ，) ，加速度a o ) ，以及前一时刻的位置f o 一剐。下一时刻的位置方程为： i o + 国) = 2 i ( f ) 一产0 一毋) + a t 2 厅0 ) ( 2 1 2 ) 上式有许多值得注意的方面，可以看到，速度并没有在式中出现。速度由f “) 的 泰勒展开式得到的附加方程解得： f o + 舀) = i 0 ) + 5 疖( r ) + 妻舀2 舀o ) + ( 2 1 3 a ) f o 一詹) = f ) 一狮( f ) + 妻国2 厅o ) ( 2 1 3 b ) 两式相减可得到速度： f o ) = 盟地杀生型 ( 2 1 1 4 ) 给定了守恒力的本算法在时间上是可逆的，保证了存储线形的动量。本方法被 证明了，即使时间步长很长也具有优良的能量守恒性质。但是，v e r l e t 方法在 速度的处理上仍然比较困难，而且此方法的形式可能不必要的会引入一些数据 的不精确。原因是在方程( 2 一1 2 ) 中，为了产生轨迹而将一个小的项( 昂2 ) 附 加到了大项( a t o ) 的微分上。 对基本v e r l e t 算法的修正解决了这个难题。其中的一种叫做半步 “l e a p f r o g 5 m 5 1 1 方法，当写下此方法后，其名字的来源就一目了然了： 第二章模型方法 机班荆+ 踟( 1 毋 ( 2 - 1 5 a ) 叭班i 卜圭国 + 筋o ) ( 2 - 1 5 b ) 首先执行速度方程( 2 一t 5 b ) ，速度通过坐标跳跃给出了下半步的值口( f + 1 2 & ) 。 在这步骤过程中，当前的速度可计算得： t o ) = ( i ( r + a t ) + t ( r 一 a t ) 2 1 6 如上式所示，l e a p 一盘o g 方法仍不能以完全满意的方式来处理速度。s w o p e ， a n d e r s e n 5 2 】等提出了“速度v e l e r t ”算法，形式为： 芦o + a t ) = f 0 ) + 西审( f ) + 妻毋2 舀( f ) ( 2 - 1 7 a ) 矿o ) = 矿( f ) + 占f 【厅o ) + a ( ，+ 古f ) 】 ( 2 1 7 b ) 首先用方程( 2 1 7 a ) 计算时刻f + a t 的新位置，半步的速度由下式计算 i ( r + 捌剐) 畦删 ( 2 - 1 8 ) 时刻h 国的力和加速度就可以计算了，速度变为 机非i ( ，+ 捌+ 圭黼( f 删 ( z 郴) 程序中速度的可视化即是直接由式( 2 - 1 7 a ) 到( 2 - 1 9 ) 的v e r l e t 方法而来的。最 后要介绍下b e e m a n 5 3 1 提出的方法： f o + 西) = i o ) + o 商( f ) + 鲁昂2 a ( f ) 一喜国2 磊o 一舀) ( 2 2 0 ) i 0 + 西) = 口o ) + 自压( r + a t ) + ；西行0 ) 一圭西面( f 一辞) ( 2 2 1 ) joo 消去速度就可以简化成方程( 2 1 2 ) 由上可知，所有的方法都有一样的全局误差，是等价的，且实际上产生了 第二章模型方法 相同的位置轨迹。 2 1 3 周期性边界条件 我们在运行计算模拟程序去预测和研究系统体积内部性质时，由于模拟的 粒子数目有限，大部分分子在样品边缘附近，因此，计算不可避免地受到表面 和有限尺寸的影响。为了减少这种影响，系统的尺寸不得不设计成极端的大， 但这样计算起来比较耗时，而且系统太大了无法模拟。因此，妥善地处理边界 条件对计算模拟非常重要。常用的边界条件有固定边界条件、周期性边界条件 等。 为了排除模拟晶胞表面的影响，在分子静力学和分子动力学模拟中，固定 边界条件得到广泛的应用。在应用固定边界条件进行模拟时，边界区域的原子 保持静止，而将可动原子局限在固定原子之内，这样就可以将模型的大小控制 在一定的范围之内。根据研究的问题不同，边界区域的原子完全固定或部分固 定。固定边界条件通常用来模拟与短程应力场( s h o r t r a n g es t r a i nf i e l d ) 有关的 点缺陷、小点缺陷团( c l u s t e r ) 的研究。但研究与长程应力场有关的问题时，如 与位错有关的问题，不得不考虑固定边界条件的应用，因为源于固定边界的虚 拟力对模拟结果产生很大的影响。 在分子动力学模拟中，当系统包含很少的原子时，原子与容器壁相互作用， 系统的行为将受到表面效应的影响，而不能产生关于犬块体积材料( b u l k m a t e r i a l s ) 的信息。因此，为了消除表面影响，合适的边界条件是非常重要的。 可通过周期性边界条件来除去这种影响。粒子被装入一个中心盒子中， 这个中心盒子通过在三个坐标方向上严格平移被无限重复。如果一个粒子在位 置f ，我们假定这个粒子实际代表一套无限个的粒子，它们的位置为 f + a + ，押云+ 芒，( ，m ，胛= 一o o ，栅) 其中z ，m ，n 是整数，西，占，手是相应的盒子的三条边在三个笛卡儿方向 ( c a r t e s i a nd i r e c t i o n s ) 的矢量。当中心盒子中的一个粒子移动时，它在盒子中 的周期镜像也以相同的方向和方式移动，因此当一个分子离开中心盒子时，它 的一个镜像将通过相对面进来，中心盒子中的粒子数目可以保持为常数。如 第二章模型方法 图2 1 所示，当一个粒子穿过边界离开时，它的所有相应的镜像也都穿过它们自 己相应的边界离开。因此在模拟中，没有必要储存所有镜像的坐标，只用储存 中心盒子中的粒子坐标。 在模拟中心盒子中，每一个粒子i 不仅与盒子中粒子，相互作用，还要和近 临盒子中的它们的镜像粒子( i m a g e ) 相互作用，这意味着相互作用能够穿过模 拟盒子边界。通过这种方式，表面影响能够被消除，既然势是短程的，周期性 边界条件的应用也不会增加计算量。 1o 一o 一o - o _o _- 0o + _ o _ o 一o 图2 1 周期性边界条件粗线组成的区域表示中心盒子 对于包含个粒子的模拟胞，用舅来表示集合粒子j 和它的镜像，那么 模拟系统中所有的粒子和它们的镜像为s s 。粒子的势和力主要来源于它 i = l 的近临粒子的相互作用，对于短程力来说，当两个粒子之间的距离大于截止距 离髓时，可认为他们将不发生相互作用。模拟盒子的三个方向上的尺度都应大 于2 尺。，这样在子集s 中任意两个粒子( i 和它的镜像) 的距离都大于截止距离； 当不同子集的粒子i 和，之间的距离小于截止距离如时，粒子i 和粒子j 的所有 镜像的距离都大于r 。，当计算作用在粒子f 上的力时，只用考虑每一个子集岛 第二章模型方法 中离它最近的粒子，而粒子，的所有镜像可以忽略。因此，利用周期性边界条 件，极大地减少了计算的粒子数。 本文的周期性边界条件采用l i n kc e l l 方法i ”】，模拟体系分成小的模拟核， 核的大小略大于原子相互作用的作用半径。在晶格中，每个原子的相邻原子有 2 6 个，在程序进行计算时仅考虑这2 6 个原子的影响，其他原子的影响通过截断 半径和e a m 势来考虑，这样就大大减少了计算量，加快了计算速度。 2 1 4 恒温的控制 在热力学平衡时，系统中每一个自由度平均都有相同的动能妄七。r ，这里 b 是玻尔兹曼常数，r 是绝对温度。系统的有效温度可以表示为动能的全体平 均( e n s e m b l ea v e r a g e ) 。 丁划去 足去 z z ， 是系统原子的数目，玎砂是瞬时温度，v i 是原子i 的速度。 在模拟中，原子的起始动能通常远离平衡组态，这使得有更多的势能转换 为动能，从而系统温度高于所要达到