Si晶体中位错运动的分子动力学模拟（可编辑）

Si晶体中位错运动的分子动力学模拟（可编辑） Si晶体中位错运动的分子动力学模拟 哈尔滨工业大学 硕士学位论文 Si晶体中位错运动的分子动力学模拟 姓名:刘艳琼 申请学位级别:硕士 专业:固体力学 指导教师:孟庆元 20051201 哈尔滨工业大学工学硕士学位论文 摘要 完全弛豫且具有低位错密度的SiGe,Si虚层衬底技术近年来在微电子， 激光，微波等领域应用的非常广泛，对其性能的研究具有重大的理论和现实 意义。分子设计和分子模拟是近来随着材料科学与计算机技术的发展而发展 起来的新研究领域，其应用前景广泛，己受到很多研究者特别是材料研究者 的重视，分子动力学模拟是分子模拟的一种常用方法。基于位错理论进行位 错动力学数值模拟研究成为当前一个崭新的研究点。本文利用分子动力学技 术对Si晶体中60度位错的运动进行了数值模拟。 本文介绍了分子动力学模拟的原理及方法，位错的运动等基本理论。根 模拟了不同温度、不同应力状态下si中60度位错的运动特性。通过使用自 行编制的分子动力学计算机程序，在Visual C++平台上开发,些代码来实 现模拟过程，并运用Matlab与VisualC++的接口技术画出位错运动速度的 变化曲线。结果明:当施加的应力大于Peierls应力时，位错才开始运 动，并且Peierls应力随着温度升高而减小;在低温时，温度的变化对位错 运动速度影响不大，而当温度升高到一定值时，位错速度随温度升高急剧减 小;当施加的剪应力小于1500MPa左右时，位错运动速度随温度的升高而 明显减小，而当应力大于1500MPa左右时，速度随应力的增大变化缓慢， 而且逐渐趋于一个极限值。温度越高，位错运动的极限速度越小。 关键词 分子动力学模拟;60度位错;位错运动速度;Peierls应力 堕堡鎏三些奎兰三兰堡圭兰堡堡兰 Abstract The with relaxed SiGe,Siheterostructureshasbeen fully technology perfect usedinthesefieldssuchas andmicrowavein widely microelectronics，laser recent theirmechanicsbehaviorsis of years(Therefore，studying great the tO academicand andsimulation meaningful its。application， Moleculardesign aNew hasa futurefor andis is field，whichbright applicationdeveloping quickly withthe ofmaterialscienceand development computertechnology(Many material been attentiontothis pursuers，especiallypursuers，havepayinggreat simulationisa methodofmolecular technique(Moleculardynamics frequent simulation(Basedonthe of simulationonits dislocation，numerical theory behaviorisan field(Inthis numerical dynamics updatestudy thesis， the simulationonthemotionofSi60。dislocationsis themolecular performedby dynamics( Inthis basic asthe ofdislocationmotion paper，sometheories，suchtheory molecular its and simulationand to dynamics methods，are introduced(According its ischosento the describeatomic characteristics，Stillinger―Weber potention interactionsin is toexertshearstresson Si;Parrinello―Rahman adopted super of under motioncharacteristics60。dislocationsdifferent cell;the temperature andstressissimulatedMolecular by process isrealizedonthe of ofmolecular platformVC十+throughcomputerprograms developedbyUS。Andthecurveofdislocationisdrawn dynamics velocity by theinterfaceofVC++andMATLAB(Theresultsindicate:Whenthe applied is than stress thePeierls dislocationto the stess，the Move，and larger begins Pererlsstressdecreasesasthe low temperatureincreases;At relatively dislocationisnot sensitivetothe temperature，thevelocity very temperature the reachsome variation，however,when temperature value， thevelocity decreasesasthe thestressislowerthan rapidly temperatureincreases;When about1500 versusthestressincreases MPa，the velocity quite rapidly，however, varies when the slowlythestressis thanabout1500MPaand velocityvery larger its the is eventuallyapproachesplateau higher velocity，The temperatureapplied， theless iSobtained( plateauvelocity (II- 哈尔滨工业大学工学硕士学位论文 molecular dislocation Keywords dynamicssimulation;60。dislocation;the Peierls stess; velocity;the -III 哈尔滨工业大学工学硕士学位论文 第1章绪论 1(1课背景及发展现状 在最近几年，基于位错理论进行位错动力学数值模拟研究成为当前一个 崭新的研究点，即微观尺度的位错动力学模拟研究和原子尺度的分子动力学 模拟研究。位错动力学模拟，属于微观尺度的研究范畴。由于它在本质上实 现了材料研究从原子尺度到宏观连续尺度的转换，包含了位错结构的演化、 第二相、以至晶粒行为研究，因此在材料塑性行为、金属硬化和材料的脆性 到韧性的转化研究中具有重要意义。这一领域的研究工作在全世界仍处于起 步阶段，主要研究小组大概有四个，美国有三个 如Zbib小组，Schwarz小 组和Ghoniem小组 ，以及法国的Kubin研究小组，国内目前这一领域的研究 工作很少见报道，周国辉【11等以A1作为研究对象，采用EAM势进行分子动 力学模拟。在I型和II型加载条件下，研究了温度和取向对发射位错的临界 应力强度因子的影响。 完全驰豫且具有低位错密度 t06,cm2 的SiGe,Si虚层 衬底技术近年来发 展极为迅速，在微电子，激光，微波等领域应用得十分广泛m】。当外延层很 薄时，外延层为了和基底的晶格结构匹配，内部就产生一定的应力。当外延 层厚达到临界厚度pJ时，在界面处就会出现失配位错来释放应力。失配位错的 出现会影响装置的性能，要提高装置的性能要尽量消除位错。但是由于si、 Ge之间的4(2,的晶格常数差，使得要得到完全没有位错的SiGe,Si结构几乎 是不可能的。通常采用配比渐变法【4】，低温硅生长法【5】或者柔性衬底技术【6】来 得到低位错密度的SiGe,Si薄膜。SiGe,Si结构的失配应力的释放主要是通过 J。 1ll 面上生成可以滑移的60度位错来实现的17 对于失配位错，在理论和实验方面已经进行大量的研究，尤其是位错的 产生和保持平衡状态时的原子结构的研究。目前，对位锚动力学行为的模拟 研究也很多。采用分子动力学研究晶体中的位错运动，对定量预测位错运动 JP等[9,10l 速度和理解晶体中复杂的塑性变形的原子机tNIsl都是有益的。Chang 了Mo中的刃位错的运动特性。线弹性动力学理论中，位错是不能超过声速 的，但是Gao等[11,121采用分子动力学模拟得到了超音速位错。另外也有人模 哈尔滨工业大学工学坝二J:学位论文 拟了金属材料铝和镍中位错的运动[1引，得到了在不同温度施加不同作用力下 铝和镍中的刃位错和螺位错的运动速度。而像Si这种晶格结构的材料中的位 错在大范围变化的温度和外力作用下的运动的模拟还未研究。 本文采用分子动力学方法，模拟在不同的温度和高应力状态下，Si中60 度位错运动的情况。研究剪应力和温度对位错运动的影响，以助于更好的理解 半导体材料的特性。 1(2分子动力学方法 1(2(1产生背景 分子模拟是一个广泛的概念，一般来说包括基于量子力学的模拟和基于 统计力学的模拟，前者为计算量子化学 computationalquantumchemistry，CQC ( 后者主要分为两个方法，分别是分子动力学模拟 molecular dynamics，MD 和 蒙特卡洛模拟 MonteCarlo，MC 。三者中以计算量子化学的结果最为可靠， 但是其计算量也是最大的，通常处理的体系也是比较小的。MC和MD都是基 于位能函数的模拟，不同之处在于MD模拟过程与时间相关，除了和MC一 样可以处理平衡性质以外，在处理传递性质等与时间相关的问题时有天然的 优势，当然MD和MC相比程序的复杂程度要高，计算的难度要大一些。 目前化学工业受关注的新技术包括单元操作集成技术、表面及界面技 术、膜技术、超临界技术、纳米技术、生物化工技术等。这些技术涉及聚合 物、两亲分子、电解质、生物活性分子等复杂物质，临赛和超临界、液晶、 超导、超微等复杂状态，界面、膜、溶液、催化等复杂现象。在这些方面除 了准确的物性数据外，更要对各种复杂现象的机理有深刻了解。 而一方面目 前的分子动力学模拟技术已经可以方便地得到许多有用的性质，例如在热力 学性质方面，可以得到状态方程、相平衡和临界常数;在热化学性质方面， 可以得到反应热和生成热、反应路径等;在光谱性质方面，可以得到偶极振 动光谱等;在力学性质方面，可以得到应力应变关系、弹性模量等:在传递 性质方面，可以得到黏度、扩散系数、热导等:在形状变化信息方面，可以 得到生物分子在晶体表面附着的位置和方式等。通过计算量子化学手段可以 准确得到物质的生成焓、偶极矩、键能、几何构型、电荷分布、各种光 谱性质等。另一方面，随着计算技术的发展，在20世纪困扰理论化学和物理 学家多时的求解问题也得到了相当程度的解决。在计算量子化学方面，Pople 喻尔滨工业大学工学硕士学位论文 等人的工作最为引人注目，在1997年获得了诺贝尔化学奖，分子模拟的软件 POLY、CHAR2MM、 也在最近的十年中有了长足进步，Cerius2、DL Gromacs等软件都已经发展得相当完善，并得到了广泛的应用。这些计算技 术成熟使人们对分子模拟的理解不仅仅停留在计算物性、补充实验手段上， 更多地可以作为一种新兴的研究手段运用在机理研究、产品设计等方面。 分子动力学 molecular 一种研究手段，目前已经有多本专著对这一学科的发展进行了总结[141。最早 对惰性气体分子，随后Harp和Berne率先提出双原子分子的模拟。对多原子 方法可以处理较小的分子。Ciccotti等提出的约束法则比较适用于更大、柔性 更强的分子体系。对于离子体系，涉及长程力的处理，最常用的处理方法是 Ewald加和、反应场方法和PPPM方法。MD模拟的系综默认为微正则系综 改进MD模拟方法。目前MD模拟本身的技术已经发展得比较完善， 研究工 作主要致力于缩短模拟时间、扩大模拟体系以及将模拟应用到不同的实际体 系中去探讨和解决问题。 1(2(2应用前景 随着分子模拟方法和高性能计算的高速发展，分子模拟已经可以运用于 很多重要的问题研究中。这种优势在于:在进行昂贵的实验合成、表征、加 工、组装和测试之前先利用计算机进行材料的设计、表征和优化，理论和模 拟可以预测出目前实验条件所无法测出的结果并可以对整个新材料的合成、 设计进行高效周全的思考。 分子动力学是计算机实验的一个分支。当今，计算机实验在科学研究中 扮演着十分重要的角色。针对量子点生长过程，分子模拟体现出以下几点优 势: 1 实验环境参数易改变。在真实的实验条件下，改变参数需要完成一系 列操作，而且需要小心谨慎。而在计算机实验中，只需要改变几个程序变量 即可。 2 观察实验细致入微，实验现象可逆。这是计算机模拟实验最突出的优 哈尔滨工业大学工学硕士学位论文 点。量子点是微小的结构，在真实实验中，观察它往往要依靠先进的显微 镜，包括原子力显微镜、移动电镜等。利用这些器材拍摄的照片，辅之以一 定的实验经验和物理功底，才能分析出结果。即使如此，仍然不能得到令人 满意的微观尺度上的图景。而采用模拟的方法，可以看到原子打入系统的实 时动画效果，生动清晰的原子表面照片，看到几十个原子整齐地排列在晶格 位置上的情况。另外，量子点的生长过程是一个不足2分钟的过程，这在人 们的感觉中非常短暂。而原子一表面撞击是一个飞秒尺度上的事件，2分钟内 X 包含了1(2 10”个飞秒，因此，这2分钟内发生的许多现象都被忽略了。虽 然已经有报道，可以用最新的实验设备拍到飞秒尺度内的实验照片。但在大 多数实验室条件下，这一点很难办到。在分子动力学模拟中，可以设定足够 小的时间步长，观察到小于1飞秒的时间间隔前后，系统的动力学行为。尤 其重要的是，在实验观测的时间范围内，我们可以操作程序，使系统向前， 向后跳转特定的时间间隔，观察系统在某一时间点上的状态。因此，由计算 机模拟得到的实验现象是可逆的，因为其本质上是程序。 3 成本较大的实验，比如高真空度、高压、极端条件下的 实验等，如果 改用分子动力学模拟的话，则会非常有益。 当然，人们会很自然的怀疑这种模拟的可信度。不过，近年来发展起来 的多体势函数，高性能计算机，新算法确保了分子动力学模拟具有足够的精 度。而且，利用分子动力学模拟作为研究手段的同时，并不意味着抛弃传统 的实验手段，把二者结合起来，分子模拟将发挥更大的作用。 1(3位错理论的研究 位错这一概念的提出是在本世纪30年代初的事情，然而由于它成功地解 释了以前所无法解决的现象，引起众多科学家的重视，从而对位错的研究提 高到很重要的地位。经科学家们的努力，创立了位错理论。今天，位错理论 己成为金属力学性能的理论基础，它是材料科学研究中不可缺少的理论基 础。位错的发展史并不漫长，但其发展过程却是迂回的。 1_3(1位错概念的提出 在早期的晶体X射线衍射强度的研究中就已发现在几乎完整的晶体中存 在有缺陷。早在1914年，达尔文 c(G Darwin 提出了图像不很明确的镶嵌组 织，用来解释实际晶体的X射线衍射强度和理想的完整晶体的差异。镶嵌组 哈尔滨工业大学工学硕士学位论文 织的概念长期为人们所沿用，但它到底是什么却无人知晓，因而它没有获得 进一步的发展。1920年，机械工程师格里费斯 Griffith 发表了一篇关于断裂 判据的文章，提出对于存在裂缝的薄板受到拉应力时，随着应变增加，裂缝 扩大，晶体强度就降低，与此同时，物理学家泰勒 GI(Taylor ]E在做铝单晶 的拉伸曲线，发现随着应变增加，晶体强度是增加的。这一结果与格里费斯 的结论恰好相反，当时无法解释这个矛盾，又一次遇到了问题。 自20年代末起，人们对金属单晶的范性形变不展了系统的研究。1926 年，弗兰克 J(Frankel 按照晶体范性形变时通过滑移面整体滑移这样的概念， 算出完整晶体的理论切变强度, 肛,s ‖为切变模量 ，而当时用实验测得 的切变强度的值仃。 10。--10。4 “，也就是说，实际晶体的屈服强度比根据 理想的完整晶体所作的理论估计值约低1000倍左右。这一现象再一次引起人 们的重视。而正在同一时期，根据理想晶体的点阵动力学理论很成功地说明 了晶体的热力学性质和弹性，这个差异更加突出，为了解释这个差异，在 1 假设，即认为晶体中存在有一种线缺陷――位错，晶体在切应力作用下内部 的位错容易滑移，并可以引起范性形变，该位错后来被称之为刃位错【”]。 1(3(2位错理论的形成 位错概念提出以后，康托洛娃 T(A(Kontopoba 又提出了一种动态的位错 点阵模型。1939年柏格斯 J(M(Burgers 又提出了螺位错的概念，把位错概念 加以普遍化，并发展了位错应力场的一般理论。接着位错理论得了多方面的 发展，并被人们用来解释各式各样的范性形变问题。从1934年位错的提出之 后，位错经历了大约15年左右的理论发展，但由于没有进行位错的直接实验 观察，当时关于晶体中是否存在位错，尚未获得实验的证明，而 对于晶体中 的位错分布情况更是一无所知，因而在一部分解释范性形变的位错理论中往 往也带有一定的任意性，从而导致了使用位错概念的金属物理学家和把位错 概念看成是“为挽救一个不成功的理论而引进一个复杂的虚构”的许多冶金家 之间的分裂，当时一部分科学家对位错理论持怀疑和非难的态度。 自1949年以后，位错理论的发展进入了一个新时期。1949年，科基尔fA H 团 ;弗兰克的螺位错促成晶体生长的理论获得令人信服的证实，多种的实验 观察揭示了晶体中位错的分布状态，证明了晶体中确实存在着位锘。位错存 在的第一个“证据”大概是从1940年的Bragg的肥皂泡实验得到的，而后许多 哈尔滨工业大学工学硕士学位论文 人几乎同时独立地在显微镜下观察到位错的存在和形状。随着透射显微镜技 术的发展，人们有可能对位错做更进一步的证实和研究。特别是在1956年以 后，人们对位错理论的了解进一步加深。1956年鲍曼 Bollmann 在不锈钢 中，赫许 PB(Hirsch 在铝中独立地发现用透射电镜透过薄到约102纳米的金 属膜能够直接观察到位错，而且由于在膜中的温度梯度或氧化层生长所引起 的应力作用下或者适当的外加应力的作用下，位错能够自由运动，Hirsch在 6J。这 80年代初已制作了位错运动的影片以及其他一些关于位错的实验结果【l 一切对于晶体中位错的结构、分布、动力学性质以及位错与范性形变的关系 等提供了确切可靠的第一手资料，证实了位错理论的一些基本论点及许多细 节，为进一步发展范性形变的位错理论奠定了巩固的基础。至此，位错理论 已经趋于成熟。 1(3(3位错理论的发展 从目前来看，位错理论的骨架已经确立，这不仅表现在位错理论已完全 能够解释范性形变中的力学问题以及位错能够通过大量实验直接观察到，重 要的还在于位错理论的确立促进了其它理论的发展，例如早在1940年柏格斯 与布莱格就提出了晶界的位错模型，在50年代这个模型得到了丰富的实验资 料的证实，从而促进了晶界理论的发展。目前，位错理论的应用范围也日益 推广，它已不仅仅局限于传统的范性形变等问题，在滞弹性、断裂、晶体的 电磁性质、晶体的光学性质以及超导体等领域内，位错理论也显得越来越重 要了。 然而，尽管位错概念是在30年代初引入晶体，位错理论知识在二战后才 获得了大的发展，历史较短，因此该理论中还存在一些不完善之处，正如弗 兰克和斯蒂兹fJwSteeds 在1975年写的一篇名为“晶体位错”的评论中所 说:“位错理论可分成不同状况的几部分，有些结论是确切的，因为它们是纯 几何的或纯形貌的，有些部分虽然是近似的，然而是非常可靠的。”例如，热 力学平衡的晶体中没有位错的这个结论，虽然对位错的热力学性质只做过粗 略的估计，但误差不会太大，不足以推翻这个结论。因此，现在有意义的问 题是不能确信那些已做的近似的可靠性，因而必须依靠全部的理论方法，观 察和推测来谋求进一步发展，除了这些“近似”之外，在位错领域中迄今还没 有完全解决的主要问题是如何填补单个位错的性质和位错集团的行为之间的 鸿沟。因此，位错理论尚有待今后进一步发展与完善。 哈尔滨工业大学工学硕士学位论文 1(4本文的研究内容及结构 本文的组织安排如下: 第1章，绪论。主要介绍了该课题的应用背景，研究意义及现阶段国内 外相近领域做出的研究成果，并提出本文所研究的内容;结合分 子动力学应 用的领域，阐述了分子动力学方法的产生、发展及其应用前景:从位错的提 出，位错理论的确立，发展演变和当前的研究领域动向上介绍了有关位错理 论研究的发展过程。 第2章，分子动力学模拟的原理及方法。针对本课题研究内容，介绍了 分子动力学模拟的基本原理和算法，包括边界条件，基本方程式，Gear预 测一矫正法，麦克斯韦一波尔兹曼分布，势函数， Parrinello(Rahmalfl算法， 并概括了分子动力学模拟解决问题的一般步骤。 第3章，si的结构及位错运动。介绍了si晶体的晶格结构， 通过论述位 错的几何形态、柏氏矢量等概念了解了位错的结构和位错运动的情况，结合本 文的研究内容，提出了建立si中60度位错运动模型的基本思路和方法。 第4章，位错运动的分予模拟及结果分析。介绍了模型的建立、模拟及结 Rahman方法旆加剪应力，用分子动力学方法模拟了不同应力和温度状态下si 中的60度位错的运动特性，通过模拟计算的结果分析位错运动的特性。 最后，对本课题的研究内容及结果做了一个全丽的总结。 哈尔滨工业大学工学硕士学位论文 第2章 分子动力学模拟的原理及方法 分子动力学方法的出发点是物理系统确定的微观描述，这个系统可以是一 个少体系统，也可以是一个多体系统，这种描述可以是哈密顿描述或拉格朗日 描述，也可以是直接用牛顿运动方程表示的描述。分子动力学方法使用运动方 程来计算系统的性质，结果得到的既有系统的静态特性，也有动态特性。 2(1边界条件 分子动力学方法中，使用有限的原子数来模拟实际晶体中原子的运动，需 要考虑表面对晶体结构中原子运动的影响。为了避免这种影响，一般采用周期 性边界条件，如图2-1所示即为三维周期性边界条件，或者固定 距离模拟原胞 表面一定厚度范围内的若干原子。在建模时，需要按照材料的晶 体结构给出其 o o o 1 o o o o o o o o o o 、o o o o,’o ? a 缀 奄 o o o o a o o o o o o o 图2(1三维周期性边界条件 原子的初始位置。为了减少分子动力学模拟在系统中粒子数小于 真实系统中粒 子数而带来的尺寸效应，分子动力学模拟采用了周期性边界条件 fPeriod BoundaryCondition，简称PBC 。所谓周期边界条件就是将一定 数量的粒子? 集中在一定的容积V中，该容积y称为元胞 SuperCell ，元胞 的基本尺寸为 上，元胞周围的部分可以看作是元胞的复制，称为镜像细胞。这些镜像细胞的 尺寸和形状与元胞完全相同，并且每个镜像细胞所包含的?个粒子是元胞中粒 哈尔滨工业大学工学硕士学位论文 子的镜像。这样，元胞在各个方向上进行周期复制便形成了宏观物质样本，因 此只需要根据元胞周围的边界条件便可以计算元胞内粒子的运动，从而可以大 减少工作量。 除了周期边界条件外，还有刚性边界条件 软边界条件 ，即限制边界粒子 的运动，减少由于边界粒子数所占总粒子数的比例太大引起的边界效应。本文 采用周期性边界条件。 2(2基本方程式 分子动力学方法是物理系统的确定性描述，它属于确定性方法。分子动力 学是通过原子间的相互作用势，求出每一个原子所受到的力，在具体的时间步 长、边界条件、初始位置和初始速度下，对由有限个粒子组成的粒子系统建立 其牛顿运动方程，用数值方法求解，得到这些粒子在相空间的运动轨迹，然后 用统计力学方法求其微观量的系统平均，得到所需要的宏观量。 7J: 分子动力学有两个基本假设【l 1 所有粒子的运动都遵循经典牛顿运动定律; 2 粒子间的相互作用满足叠加原理。 因此，分子动力学虽然是从原子层次研究问题，但是其本身忽略了效应的 影响，这里的原子点实际上属于宏观粒子点的范畴，只是粒子问的作用特性通 过势函数体现。本质上，可以简单的将分子动力学看作是广义牛顿运动方程的 数值积分，仍然属于近似计算。 对于粒子数为tl的物理系统，其力学描述的哈密顿形式可以写为: aH 口 ―― 1 ? i l，2，3„n 2-1 P:一塑 2―2 oq， 哈密顿函数?为 ?:y旦+u ‘?2州 2-3 式中q，――第i个粒子的广义坐标 晗尔滨T业大学丁学硕士学位论 文 P，――第i个粒子的广义动量; 聃――第i个粒子的质量; U――系统总势能。 浚系统力学描述的牛顿形式则为 l f?J,m， i-1，2，3„n 2-4 E 一V。u“，r2，r3，„oJ 2-5 式中一――第f个粒子所受的合力; n――第i个粒子的位置坐标; [，――势能函数。 经典力学方程还有其它一些形式，虽然这些方程在数学 上是等价的，但在 数值上却是不同的。为简便起见，常选择牛顿方程形式。 2(3 Gear预测一校正法 牛顿运动方程的一般形式为: m，l f 2―6 式中m(――单个原子的质量: F――由式 2-5 得出。 通常我们求解方程 2-6 的方法是用有限差分法，即己知原子在r时刻的位 置、速度和其他动力学信息，求t+6t时刻原子的位置、速度等量，并使它们 满足一定的精度要求。毋的选取依赖于特定的算法，而且必 须远远小于分子 移动一个与它自身尺寸相当的距离所需的时间。有许多算法可以用来解决这样 的问题。如果经典粒子的轨迹是连续性的，我们可以用Taylor展开的方式得到 ，+8t时刻粒子的位置和速度等信息。 Prt6++ p+毋 ，o + f+毋 哪 + +去5r36p +-一(( 6tv t 6t2a t +扣? +????? aP t+6t、 a t +6tb t +(((( bP t+6t 6 ， +„( 别代表一组坐标位置和速度值。a是加速度，b是，的3阶导 数。如果我们裁 哈尔滨工业大学工学硕士学位论文 断Taylor展开式，就把它写到如式 2(7 所表示的那样。那么从表面看来，我们 从当前粒子信息推断出了下一个时刻粒子的信息。在这样的做法中，我们需要 把，，n“，b存储到计算机中。其它的算法也许会要求计算机记住不同的 加速度a t一& 等，并得到不同于式 2―7 的预测方程。但是，不管怎样，方程 f2―7 不能得到足够精确的结果。因为到目前为止，我们还没有把牛顿运动方程 引入进来。所以我们由预测值，’ t+fit 求出t+fit时刻的力和校正加速度 “。f，+占， 。这样得到的加速度再和预测值相比较，产生 一个误差信号: Aa t+西 ac f+研 一日’ f+西 2―8 这个误差信号和预测值一起带入校正阶段的运算中，通常的表达式如下: rC t+fit rP t+6t +coAa t+fit 2_9 8t + mvc t+fi?t 纠ve一 t+帅St +t口‘ 十Jr 口’U+ + clc2幽Aao tZ8三t bc t+fit be“+60+c3Aa t+dt 这样能够得到对真实速度、坐标等更好的近似。Gear曾经讨论过系数?, o等的最佳选取问题。如果预测或校正的变量不同，运动方程对r的阶数不 同，或是同阶微分方程的具体形式不同，系数的选取会有所不同。如式 2。7 、 2(9 所表示的是预测校正方法的最普遍形式，被称为Nordsieck表达方式。人 们发现使用被时间步长标定的速度、加速度等量来描述预测校正方法是非常方 便的。一般地，我们把,定义为一组粒子的坐标。并以此为参照，定义了一系 列被时间步长标定的导数值; 1 ( 1 3 2ro,dt 3，o,dt 2―10 _ fit dro,dt ，如 ?函2 d2 ，r3 ?研3 d ‘ o 举例来说，有4个标定量的Gear算法【l引，它的预测式的矩 阵形式为: l 1 1 l ，o r O l 2 3 rl r f2―111 0 0 l 3 ，2 r 0 O O l ? f 这是一个上三角矩阵，给计算预测方程带来很多好处。相应地，预测一校正方 法的校正式可写成: 哈尔滨工业大学工学硕士学位论文 r， t+研 CO ，。’ r+& C1 J- 2―12 ，2，J f+出 C2 一’0+研 C3 2阶移动方程的Gear校正系数列在表2-1中 移动方程的形式为: r f r、 2-13 如果方程的形式为 r -， ，，r 2―14 一阶导数项出现在等式的右边。则在有5个预测校正量的方法中，c。应改为 19,90;在有6个量的方法中，cn应改为3,16。本文所采用的是2阶移动方 程，形式为式f2(12 ，具有六个校正系数的Gear预测校正方法。 综上所述，Gear预测校正方法共分3步: 1 预测。由t时刻粒子的坐标、速度等量通过Taylor展开的方式得出 ，+6t时刻的值。当然，在所有这些被预测的量当中，必须包含加速度口。 2 估计力。由预测得到的坐标值，对势函数取梯度，计算力，进一步求 加速度，这样得到的加速度通常会与预测值不同，这种不同以“误差信号”的形 式表示出来。 3 校正。误差信号可用来校正位移与它的导数项。所有校正项都与误差 信号成jF比，如式 2―9 。如上文所讲述的那样，确定校正系数显得十分重要。 表2-12阶移动方程的Gear校正系数 系数总数 C0 q C2 岛 C4 C5 3 O 1 1 4 l,6 5,6 1 113 5 19,120 3,4 1 112 1112 1 6 3,20 251,360 11,18 1,6 1,60 哈尔滨工业大学工学硕士学位论文 2(4麦克斯韦一波尔兹曼分布 在进行分子动力学模拟之前，我们首先假设粒子服从麦克 斯韦一玻尔兹曼 分布 M―B分布 ”9】。M(B分布是近独立经典粒子按能量的最 概然分布。其特点 是:粒子彼此可以区分，每个量子态中的粒子数不受限制。下面让我们具体讨 论一下M―B分布。 首先，什么是经典粒子?从统计的角度看，用经典手段描述的经典粒子和 用波粒二象性描述的量子粒子的区别在于:经典粒子的状态可以用确定的坐标 和动量来描述，既它具有自己的“个性”，从而是可以区分的;而量子粒子不能 同时具有确定的坐标和向量，其运动状态只能用全同波函数来描述，通常是不 可区分的，量子粒子又可以分为费米子和玻色子。 在大量粒子组成的系统中，粒子如何按能量分布，是统计物理学的,个基 本课题。设体系由N个近独立粒子组成，体系中能量为邑，岛„(s(„(的能级 上分别有?，，?:„???，„??个粒子，每个能级包含的量子数为g。，g:，„??gf'„，体 系的总粒子数?和总能量E守恒，即 fyM N:恒量 黔班 龃 。‘15’ 我们在上述条件下计算粒子按能量的最概然分布。 首先说明等概率假设。对于处在平衡态的孤立系，其宏观态完全确定，但 系统的微观态并不确定，即不能认为系统的某一可能的微观状态比另一微观状 态出现的几率大。1871年玻尔兹曼提出著名的等概率假设，具 体表述为:对于 处在平衡态的孤立体系，其各个可能的微观态出现的几率相等，也就是说:如 果平衡态孤立系的可能微观态的总数为职则系统的任意微观态出现的概率为 I,W(即 P】 r P2 r 一’P。 r 素 2-16 等概率假设是平衡态统计理论的基础，其准确性己被大量实验所证实。根 据等概率假设，处于平衡态的孤立系，每一个可能的微观态出现的几率都相 等，在各种可能的分布中哪一种分布对应的微观状态的数目多，此分布的概率 就大。我们把出现概率最大的分布叫做体系的最概然分布。 在热力学中，一个孤立体系最终要达到热力学平衡态的，从统计物理学的 角度看，这就是系统自发的趋向于最概然分布，具体思路如下: 哈尔滨工业大学工学硕士学位论文 1 求将N个粒子按?l，N:，„??，N一-??这种分布放到能量s，，s:，„f，„的各 种量子态中去的可能占据方式数Q。 2 求Q取最大值的分布，即最概然分布。 2(5势函数 在分子动力学模拟计算中，势函数的选择非常的重要，它决定了计算工作 量和计算模拟与真实系统的近似程度。从能量的角度看，原子结 合成晶体时， 晶体的总能量应当比原子在自由状态下的总能量低，晶体才能处于稳定状态。 晶体的总能量与自由原子的总能量之差定义为晶体的结合能，即 , U一, 2―17 式中 乙L一是绝对零度时晶体总能量: u。――为组成晶体的自由原子的总能量。 很明显如自由原子能形成稳定的晶体，结合能应当是负的，其数值I,J就 是把晶体分离为自由原子所需要的能量。结合能大的晶体，原子间的结合比较 强。但是模拟结果的准确与否的关键在于材料系统内的原子间相互作用势函数 [2叫的选取。 原子涮的作用势依赖于原子间的距离和所带电荷。固体材料的诸多性质最 终都决定于原子间的相互作用，所以在固体材料的原子模拟中，作用势的选取 是计算机模拟中的一个极为重要的环节。 合理选择势函数的要求一般包括为:能够足够精确，在实验范围内，计算 结果能较好的和实验现象吻合，在实验范围外，能够提供合理的预测结果;能 够精良简单，在有保证的计算精度的前提下，应使计算机处理起来简便易行， 速度足够地快。 对于Si和Ge等半导体，其键和强度依赖周围原子的配置，必须引入三体 稳定的势函数，一般只用与描述Si和Ge等半导体材料。本文是模拟si中的 位错运动特性，因此选用S(w势函数来描述si中原子问的相互作用。 哈尔滨工业大学工学硕士学位论文 ? ?V: f，， +?V3 f，，，k ; 当 :缆t V: o s^ o,o 2一18 v3 I，，o，珞 t艿 1,盯，，‘,仃，‘,盯 式中,，a分别表示与键的强度和原子大小相关联的参数，正和 五两个函数的具 体形式如下: 五 ， :f一 矿V9 e坤[p 2-19 10，r?口 z―z， 其中，向 。，k，, Aexp[r 。一口 ,+， k一。 “]x c。s吒+; 1 式中的其他参数都与材料本身有关(确定研究材料后都可确定，本文是研究Si 晶体，由Si晶体的材料性质确定这些参数的信。 2(6 Parrinello(Rahman算法 最初的分子动力学方法只适用于模拟单元的体积和形状不变的模型中，而 我们研究的模型很多都需要在一定的温度下，有一定外力作用下的情况，这就 必然会使模型发生体积和形状上的变化。因此，将模拟单元的体积和形状维持 为常量的约束受到解决晶体结构变形问题的方法的可用性的限制，在这类变形 服了这个难题，本文在模拟si晶体中60度位错的运动中，就采用了此方法施 加剪应力。 水压力下的粒子体系在体积和形状上的变化。而在Andersen方法中，仅容许 分子动力学fMD 单元在体积上的改变，因而由于基本能量流受到抑制，即MD 单元形状受约束，晶体结构变形在Andersen方法中是受限的。 以下是分子动力学中多余弹性的引入:如前，一个可以被周期性复制填满 整个空间的单元包含?个粒子，而这种单元却可以有任意体积和形状完全是由 哈尔滨工业大学工学硕士学位论文 组成MD单元边缘的a、b、C三个矢量所刻画出的。三个矢量可以有不同的长 度和任意共同的起始点，这种描述的可变性是通过矢量 a，b，c 形成一个3×3 的矩阵来实现的，h的各列依次为a、b、C三个分量，体积由式 2(22 给出: n llhlI a?? 6，、c 2-22 a、b、c依次构成右手系，粒子f的位置坐标n可由h和列向量s。得出，带有 参数喜，叩，和,， f2向s，2善a+叩，6+,c 2-23 显见变量范围为。茎专，珥，,sl i 1，(，N。sj的模为s+ ^，,2，y ，其中五，‖和 u为 一。。，+o。 的积分。以通常的方式标量化一个矢量或张量，而后粒子 i，J之I'BJ,E离的平方表示为: 哼2【I一? G 墨一‘ 2-24 其中张量矩阵G为， G_龇 2(25 为完整说明用于此处的批注，最终我们以下面矢量度量此交互的空间: 等 6AC,CAQ，GA6 ;罟仃 2-26 矩阵盯;Qh二1包含了关于MD单元大小和原点的全部信息。 1(仅受静水压力的情况 如何获得MD单元的形状和大小的改变:通常设出3n个动态变量描述n 到: 妒 1,2?硝啦一??庐 o +1,2珊“一pQ 2―27 其中P是我们选定要加载静水雎为的系统，矽是以后要进一步阐明的一个 质量维数，是否可由第一性原理导出此Largragian式是下一步 研究的一个方 向;现在，通过运动方程和由其产生的假设可验证其正确性。 从方程 2―27 得 到: j 一?叫 妒’70 *一t 一G。峦(，i 1„N 2― 28 ‖矗 万一p 盯 r2―29 其中，用通用的相关符号，写作v 艇 哈尔滨工业大学工学硕士学位论文 Q万 ?鸭V，_-EE o’70ho r2-30 f J J l 式， J?，1”,? 当然体系中的压力并不是可控的，其值可由通常的方式取l,3的n的平均轨迹 获得。 Hooveretal导出的方程【矧与方程 2(28 类似，在应变率张量转换为ll 1。 时，它们的两个一阶运动方程是等价于方程 2―28 的，因此h- 应变率张 量 ’h，而就如所希望的，h受应变力张量的驱动，它们的方程适用于特殊的 驱动不平衡现象的研究。 方程 2―29 允许系统受到一种动态平衡的驱动，这种动态平衡介于外加应 力和内应力张量之间，因而方程 2-29 允许应用于以上提到所有平衡驱动的非 平衡现象的研究。在某种平衡态，使外应力有个振荡的时间独立也将容许用于 系统对多种外部激励的独立频响的研究。从方程 2(29 ，相应的质量‖决定了 从介于外部压力和内应力间的平衡态恢复过来的时间，正如Andersen谈到 的，对‖值的恰当选取使量级为相同阶数的驰豫的时间成为更大试样的--'b 部分驰豫，他提出一种选择是与L,c有相同量级的阶数，其中上是MD单元的 大小，C是声速，这就明显限制了矿的选择度，并且使得估算更 加接近于实 际。然而，如果只关注于静态平衡，缈就可以根据计算的便捷来选取。实际 上，在经典静力学理论中，一个系统的平衡比率是独立于其中的质量组分的。 从方程 2(27 可以推出相应的服从通常力学法则的哈密顿量。由于系统的 外力与时间不相关，这就成了一个运动常数。就有: 2??32 A ?l,2m,v2,+ZZ? r, +1,2WTr,t’五+pQ J J j i 温度为r的平衡态下，9,2br由带有‖的项和其他动态项得到的3N,2如r构 成。因此当达到3:N的精度时发现运动常数h就是焓 日 E+pf 2-33 其中 E ?1,2m,v2,+??? o 2-34 j ，』 』 因此由式 2―27 中的Lagrangian项产生一个0，E? 分量。 2(受一般应力的情况 哈尔滨工业大学工学硕士学位论文 以上简述的分子动力学公式，自然地被放到各向同性的外部 应力情况下进 行讨论。不难看出由于是在经典弹性理论中，应变的概念是与矩 阵张量形式的 变量密切相关的并且就如所显现的，G是分子动力学中一个必不 可少的部分。 为了使以上的讨论具有可行性，像弹性力学通常所做的那 样，需要引入一 个参照态。应用引入的概念，系统的参照态可由其矩阵ho和体 积Q。 lfh。||定 义。在此参考系中由矢量S给定的空间中的一点的位置 ro hoS 2-35 系统的,个均匀扰动使ho变为h(，从ro移至r，其中 r5hs2hh;。ro 2-36 根据扰动给定位移u: 萨r-Ij 1111;-1 岛 r2―37 由Landau和Lifshitz定义应变张量占，用扎表示r0分量 ? Vz 鼍+善+?薏薏] cz。s， 我们发现，根据定义G的式 2-25 和定义U的式 2(37 ，有 占 1,2 h;1G列一1 2―39 确定应变占，弹性能,的表达式就可写出。假设S是外部应力， P是静水压 力，则: 屹 p f,-f20 +QoTr s―p s 2―40 在微应变的限制下，Tre Af ,f20 Q一? ,ao 2―41 因此当式 2(41 作了正确的近似时，可以得到更通俗的表达， , QoTrSe 2-42 否则，当式 2―41 不精确时，需要通过式 2-40 雕J静水压力效应下的正确描述。 为简化式 2―27 刨J pD，这就给出了新的Lagrangian项吼， 织 妒一1,2TrEG( 2-43 其中对称张量?与应力的关系为: E h。。1 S-p h'。。Qo 2(44 哈尔滨工业大学工学硕士学位论文 等式Tr AB Tr BA 。 2―29 变为 嘲 万一p 盯一hE 2??45 中的焓的一般表达为: 凤 E+K， 2-46 其中E由式 2-34 给出，圪由式 2―40 给出。 运动方程 2(45 表示一个平衡态需要使其右端的平均值为 零。由?的定义， 从方程 2(26 的平衡写出GO，由此得到 万一p 盯 h h; s-p 盯o 2―47 这表示，就如直观上看到的，对于一个平衡系统，将常量矩阵? 与外加应力矩 阵S联立，对参照态的最合理选择取决于ho h 。 2(7分子动力学模拟的技术路线 下面让我们总结一下用分子动力学解决一般问题的做法: 1(初始化。输入原子团信息，建立势函数。原子团信息包括原子类型， 质量，晶格特征 空间点阵结构、晶格常数等 。另外，还要声明是否使用周 期性边界条件，给定有限差分法的步长，规定模拟过程中的外部条件。 2(动驰豫。在程序开始时，原子通常被固定在晶格位置上，相应于给定 的势函数，这种结构在T 0时是平衡的。因为系统的对称性，在一个建构良好 的晶格中，每个原子所受的合力都为0，原子始终处于各自的晶 格位置上，移 动方程将不能求解出粒子的运动轨迹。因此，在起始阶段，粒子应被赋予随机 速度。我们给粒子赋初值，让它们等于0或遵从M―B分布。当然，这种初始 状态并不满足平衡条件。不过，一旦程序运行起来，系统会在100个计算步长 内达到iF衡。随后用有限差分法求出粒子在每个时间步长内的位移、速度和加 速度等。有两种方法可以用来引入随机分量: f1 给晶格位置一点随机位移，位移的幅度不可过大，以免原子核重叠。 f2 赋予粒子初速度，并使它们满足特定温度下的M―B分布。当我们这样 做时，系统整体上将有一个较小的线性动量，体现出此时系统的平动。因为这 哈尔滨工业大学工学硕士学位论文 种情况是不大好处理的，通常要从每一个粒子中减去这一分量，以使系统的总 动量为0。值得一提的是，在起始阶段引入随机量是分子动力学模拟中唯一体 现随机性的地方，随后展开的计算过程都是确定性的。 3(统计物理的方法求宏观特性。统计物理学是研究热现象的微观理论。 它认为宏观现象是大量粒子运动的集体表现，宏观量是相应微观量的统计平均 值。并在对物质的微观结构、微观粒子间相互作用做出某种假设的基础上，去 分析微观粒子的力学运动，从而揭示热现象的微观实质。统计物理【25】研究的对 象是大量粒子组成的微观体系。我们把描述体系整体特征的物理量称为宏观 量。宏观量即物理量的宏观测量值。如体系的总质量，气体的压强，温度等。 而单个粒子所具有的物理量，如每个粒子的质量，瞬时速度，动量，量子数 等，我们称之为微观量。我们把不受外界影响时，系统的宏观量具有稳定值的 状态称为系统的平衡状态。这里的不受外界影响指的是系统与外界无能量交换 不传热，不做功 ，同时系统内也没有化学反应、热核反应等能量交换。当 然，绝对排除外界影响时做不到的 例如体系总是受到引力场的 作用 ，严格的 平衡状态只是一种理想情况。但即使在这种情况下也只是系统的宏观性质不发 生变化，体系内粒子的热运动一刻也没有停止。即体系的微观性质是在不断变 化的，体系的平衡状态只能是一种热动平衡。所以，统计物理认为，任何宏观 量都只能反映一系列随时间变化着的微观状态所产生的平均效果，系统的宏观 量都是在测量时间内，系统所有的微观状态中相应的微观量的统计平均值。最 后需要指出的是，在许多分子动力学的研究工作中，求统计平均值这一环节都 是必不可少的。在本文研究过程中需要确定位错芯的位置，就是通过寻找滑移 面上的高能原子，求其平均位置，再求得均方差值下来，从而确定位错芯 位置。 2(8本章小结 本章主要内容: 1 针对本课题研究内容，介绍了分子动力学模拟的基本原理和算法，包 括边界条件，基本方程式，Gear预测一矫正法，麦克斯韦一波尔兹曼分布，势 2 概括了用分子动力学方法解决问题的一般步骤:初始化，动弛豫，用 统计物理的方法求宏观特性。 堕查堡三些查耋三兰堡圭耋堡竺兰 第3章Si的结构及位错运动 3(1 Si的结构 si是我们这一问题的研究对象，所以有必要搞清它的晶格结构。Si是一种 钻石结构的晶体‘261。钻石结构的空间点阵是面心立方。它的一个基本单元由两 0 1,41,4位置上，这样的单元放到每个 个相同的原子组成，分别位于0 O;1,4 点阵点上，就形成了钻石结构。另一方面，钻石结构可以认为由一个面心立方 结构沿体对角线平移而得到，如图3(1所示。图3-2中显示出原子之间的结 图3-1钻石结构晶体中原子的相对位置 合键，一个si原子和周围4个si原子相互作用。C，si，Ge都可以形成钻石结 0 1,41,4上的 构。为了在程序中得到这种钻石结构，我们把最初位于0 O;1,4 两个si原子反复复制，具体做法如下: fill 2 particle Si 00O Si O|250(25O(25 fillcell O O(5 O(5 O(5 O 0(5 O(5 O(5 0 O(5 0 0(5 其中， 0 0(5 0(5 O(5 0(5 0 表示Si原子单元的3个移动方 向。符合钻石结构的原子投影在立方晶胞某一表面上的情况。分数代表以面心 哈尔滨工业大学工学硕士学位论文 立方晶胞的边长为单位时，原子距离底面的高度。表示为0，1,2的点位于面 心立方的点阵点上;表示为1,4，3,4的点位于另一个面心立方晶胞的点阵点 上，这个晶胞可以通过把目前的晶胞沿对角线方向移动对角线长度的1,4来得 到。 图3-2钻石晶格结构 3(2位错的结构 位错是晶体中存在的非常重要的晶体缺陷，属于线缺陷。在它周围的严重 畸变区域的尺度中，其中一个尺度上比垂直于此线的其它两维尺度大得多。也 可以把严重畸变区域用类似一个“管道”来描述，这个管道的直径通常仅有几个 原子间距，并贯穿于晶体之中。在管道内，原子间的坐标与在完整晶体中很不 同，而在管道之外的原子的坐标接近于完整晶体。这里的所谓管道“内部”和管 道“外部”这个定义不很精确的区域成为位错核心。 3(2(1位错的几何形态 任意形状的位错环的结构是复杂的，但是不论形状如何复杂的位错线都可 以看成是由两种简单类型的位错即刃位错和螺位错混合而成。在此，以简单立 方晶体为例来看这两种简单位错的几何形态口”。 1(刃型位错 图3(3表示一块单晶体，在切应力r的作用下发生局部滑移。其中ABCD 为滑移面，ADEF为已滑移区，在此面上部的晶体相对于下部向右滑移一个原 予间距，EFCB为未滑移区。发生局部滑移后晶体内出现了一个多余半原子面 EFGH，显然在晶格内产生了缺陷，这就是位错。这种位错在晶体中有一个刀 哈尔滨工业大学工学硕士学位论文 刃状的多余半原子面，所以称为刃型位错。 图3-3刃型位错 盯线是多余半原子面和滑移面的交线，也是已滑移区和未滑移区的交 界，这条线就是位错线，它与滑移方向垂直。位错在晶体中引起的畸变在位错 线处最大，离位错线越远晶格畸变越小。原子严重错排的区域约3-4个原子间 距，因此，位错是沿位错线为中心的一个管道。图3-4是仅取三个原子层厚度 的简单立方晶体中一个刃位错的原子排列模型，通常把多余半原子面在滑移面 以上的刃型位错称为正刃型位错;反之，把多余半原子面在滑移面以下的刃型 位错称为负刃型位错。 幽3-4刃型位错原子排列模型 2(螺型位错 晶体以图3―5所示方式发生局部滑移一个原子间距，ABCD面为滑移 面，EF线以右为已滑移区:以左为未滑移区，EF线为位错线，它和滑移方 向平行，沿位错线原子面呈螺旋形，每绕轴一周，原子面上升一 个原子间 距，这种位错称螺型位错。螺位错分为左旋和右旋，根据螺旋面 旋转方向， 篁堡堡三些奎耋三耋竺圭耋堡篁 兰 符合右手法则 即以右手拇指代表螺旋面前进方向