分子动力学模拟位错和界面的相互作用

1、学校代码 10530 学 号 200910081121 分 类 号 密 级 硕 士 学 位 论 文分子动力学模拟位错和界面的相互作用 学 位 申 请 人 周银库 指 导 教 师 陈尚达 副教授 学 院 名 称 材料与光电物理学院 学 科 专 业 材料科学与工程 研 究 方 向 金属薄膜的力学性能 二零一二年五月Molecular dynamics simulations of interaction between dislocations and interfacesCandidate Yinku Zhou Supervisor Shangda Chen (Associate Profess

2、or) College Faculty of Materials, Optoelectronics and Physics Program Material Science and Engineering Specialization Mechanical Properties of metal film Degree Engineering Master University Xiangtan University Date May, 2012 湘潭大学学位论文原创性声明本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含

3、任何其他个人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承担。作者签名：日期： 年 月 日学位论文版权使用授权书本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授权湘潭大学可以将本学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。涉密论文按学校规定处理。作者签名：日期： 年 月 日导师签名：日期： 年 月 日摘 要纳米尺度的金属多层膜在屈服应力、塑性、抗腐蚀性能

4、等方面具有特殊的性能。目前它已被广泛应用于航空航天、机械制造、电子技术、光学工程及计算机工程等各个领域。而在薄膜材料的应用过程中，薄膜的使用寿命和可靠性是人们普遍关注的焦点问题。界面的结合性能是影响多层膜寿命和可靠性的关键指标，而位错和界面的相互作用机理决定着界面的结合性能，即位错和界面的相互作用机理在薄膜的使用寿命和可靠性方面扮演着关键角色。因此对位错和界面的相互作用机理的研究就显得特别有价值和意义。随着高性能计算机的发展，原子模拟已成为材料性能预测与设计方面一种有效的方法。 本文用三维分子动力学方法研究了位错和界面的相互作用机理，具体如下：首先，用分子动力学方法研究了侧向拉伸载荷下位错从b

5、cc-Fe/Ni界面的形核和发射过程。弛豫后，在Fe(0 0 1)/Ni(0 0 1)和Fe(0 0 1)/Ni(1 1 1)界面观察到无序的失配位错网络，Fe(0 0 1)/Ni(1 1 0)界面观察到长方形的失配位错网络。研究了晶体取向对Fe/Ni双层膜拉伸性能的影响。不同取向的对比发现Fe(0 0 1)/Ni(1 1 0)系统的屈服强度最低。和Fe薄膜进行了对比，发现Fe/Ni双层膜系统的塑性高于Fe薄膜的，而屈服强度低于Fe薄膜的。模拟结果显示，界面是位错的发射源，滑移位错从界面的失配位错线形核和发射。同时界面也会阻碍位错运动，随着拉伸的进行，Fe层中越来越多的位错被塞积在界面处，当到

6、达到临界值时，迫使位错穿过Fe/Ni界面，从Fe层到Ni层。在Fe基体中位错主要在1 0 1面滑移，而在Ni中主要在1 1 1面滑移。其次，用分子动力学模拟了单轴拉伸载荷下不同扭转角的Cu(001)/Ni(001)界面的结合性能。模拟结果显示，当扭转角小于15.124度时，界面形成方格状的失配位错网络，界面失配位错网络的密度随着扭转角的增加而增加。当扭转角大于15.124度时，在界面形成面缺陷。模拟发现界面构型对Cu/Ni系统的界面强度有着非常显著的影响。随着扭转角的增加屈服应力首先减小，直到扭转角为5.906度的最小值，然后增加，当其达到扭转角为15.124度的最大值后，又开始减小 ，最后当

7、扭转角约大于20度，屈服应力几乎趋于一稳定的值。关键词： 分子动力学；界面；滑移位错；失配位错 ABSTRACTNanoscale multilayered composites often possess extraordinary mechanical properties in terms of yield stress, ductility, and wear resistant. Now it has been widely used in aerospace, mechanical manufacturing, electronics, optical engineering an

8、d computer engineering fields. In the applied process of the thin film materials, people commonly focus on the reliability and service life of the thin films. The binding property between the thin film and the substrate is a key indicator of metallic multilayers the reliability and service life, the

9、 interaction mechanisms between dislocations and interfaces dominate the binding properties of the interfaces. Therefore, the interaction mechanism between dislocations and interfaces plays a vital role in the field of the reliability and service life of the thin films.So it is an interesting and va

10、luable thing to understand the interaction mechanisms between dislocations and interfaces. With the development of high-performance computer, atomic simulations have become an effective method in the field of material properties forecast and design. In the present work, we have studied the interacti

11、on between dislocations and interfaces with 3D Molecular Dynamic Simulations. Firstly, molecular dynamics simulations were carried out to investigate the nucleation and emission of dislocations from an interface in a bcc-Fe/Ni bilayer subjected to transverse loading. After relaxation, disordered typ

12、es of dislocations were observed at both Fe(0 0 1)/Ni(0 0 1) and Fe(0 0 1)/Ni(1 1 1) interfaces, and rectangular dislocations types at Fe(0 0 1)/Ni(1 1 0) interface. The orientation effect on the mechanical properties of a Fe/Ni bilayer system was investigated. The yield stress of the Fe(0 0 1)/Ni(1

13、 1 0) system abtained is lowest. We also found that the yield stress of pure iron nanofilm was higher than that of a Fe/Ni bilayer system, and the ductility was lower than that of a Fe/Ni bilayer system for given temperature and strain rate. The simulation results obtained also show that the misfit

14、dislocations at Fe/Ni interface acted as a source to nucleation and emission of glide dislocations. Glide dislocations nucleation and emission from misfit dislocation line at Fe/Ni interface. The existence of misfit dislocations and the lattice mismatch can also act as barriers to dislocation motion

15、 and transmission across the interface. More dislocations in Fe have been arrested at the Fe/Ni interface, which provides sufficient stress for dislocations to transmit from Fe to Ni. Glide dislocations mainly occurred on 1 0 1 plane in Fe layer of FeNi bilayer, and 1 1 1 plane in Ni layer. Secondly

16、, molecular dynamics simulations were carried out to study the mechanical properties of Cu(001)/Ni(001) interface boundaries with different twist angles subjected to uniaxial loading. The results obtained revealed that square misfit dislocations networks can be observed when the twist angle was lowe

17、r than 15.124°, and the density of misfit dislocations increased with increasing twist angle. Face defects were formed when the twist angle was higher than 15.124°. It has been found that the interface configuration had a significant effect on the interface strength of the Cu/Ni system. Th

18、e yield stress was found to decrease first with increasing twist angle and it reached its lowest value at 5.906° twist angle. Subsequently, it increased with increasing twist angle till it reached its highest value at 15.124° of the latter; it then decreased again and finally became almost

19、 constant when the twist angle was larger than approximately 20° .Key Words: Molecular Dynamics; Interface; Glide dislocation; Misfit dislocation 目 录第1章 引言11.1 薄膜概述11.2 界面概述11.3晶体位错相关理论概述41.4晶体界面的分子动力学研究现状61.5本文的研究思路、目的及意义7第2章 分子动力学方法92.1 引言92.2基本原理102.2.1 积分方法102.2.2 原子间的相互作用势112.3 边界条件132.4 温

20、度、压力控制方法142.4.1 控温方法142.4.2 控压方法152.5 分析方法162.5.1 中心对称参数（centrosymmetry parameter）162.5.2 径向分布函数172.6 模拟软件与可视化172.7 小结18第3章 Fe/Ni双层膜系统中位错和界面的相互作用193.1 模拟方法193.2 结果及讨论203.2.1 弛豫后沿厚度方向的势能分布203.2.2初始界面失配位错结构213.2.3 晶体取向的影响223.2.4 bcc-Fe薄膜的变形机理233.2.5 bcc-Fe层的滑移位错243.2.6 Fe/Ni双层膜的变形机理263.2.7 Ni层厚度的影响293

21、.2.8 温度的影响303.2.9 应变率的影响303.3 小结31第4章 Cu(001)/Ni(001)扭转界面的结合强度334.1 模拟模型和方法334.2 结果与讨论354.2.1 初始失配位错构型354.2.2 Cu(001)/Ni(001)系统的变形机理364.2.3 Cu层中滑移位错的变形机理384.3 小结39第5章 总结与展望415.1 工作总结415.2 工作展望41参考文献43致 谢48个人简历、攻读硕士学位期间发表的论文49V湘潭大学硕士毕业论文第1章 引言1.1 薄膜概述薄膜材料通常是指在二维方向上的尺度比另一维方向大很多，在其基体材料上通过化学或物理等方法制备另外一层

22、材料，以达到某种特定功能，如提高断裂韧性、增加塑性、增加使用寿命, 以及美观度等等。薄膜材料与块体材料相比，有比较高的比表面积，比较少的配位数，由于原子之间的化学键在表面处突然断开，故使其表面原子相对内部原子有较高的表面势能和表面应力1-3。有很多种薄膜的制备方法4，如物理方法中的真空蒸发、溅射、离子束和离子助、外延膜沉积技术等，以及化学方法中的电镀、阳极反应沉积法、化学气相沉积等。各种制备方法有各自的优缺点和适用范围，实际应用时以具体情况而定。薄膜的分类方式有多种5。按性质来分，有合成膜与天然膜等; 按层数,有单层膜和多层膜等。由于具有独特的微观结构以及物理和化学性能，薄膜材料在最近几十年来

23、已被广泛应用于航空航天、机械制造、电子技术、光学工程及计算机工程等各个领域。目前许多国家都把薄膜材料的研究设为大型研究项目。在各种类型薄膜中，金属多层膜在屈服应力、塑性、抗腐蚀性能等方面具有特殊的性能6-11。既而金属/金属薄膜材料被广泛应用于各种结构功能材料中，薄膜材料应用的可靠性和使用寿命很大程度上依赖于薄膜与基体的结合性能，据统计，日、美、欧共体等国每年国民生产总值的6%8%12 都因材料的疲劳、断裂、腐蚀、磨损等破坏而损失。而界面的结合性能是影响多层膜质量的关键指标，这就要求我们对薄膜界面微观结构有很好的了解。对界面结合性能的研究，已经是当今具有困惑性和挑战性的难点问题。目前已有大量的

24、实验研究方法，如纳米压痕法、鼓包法、拉伸法等。另外，随着近年来高性能计算机的发展，计算机模拟方法也被越来越广泛的应用于这一领域。1.2 界面概述界面一般包含表面、晶界和相界面等，为一种二维缺陷，即面缺陷。我们称晶体与空气或液体接触的界面为表面。在表面上的原子，其相邻原子数比晶体内部要少，相当于一部分结合键被拆断，因而有较高的能量，产生了表面能，表面能通常要比晶界能更大些。通常应用中的金属材料，大部分为多晶材料，而非单晶材料。在多晶材料中，在晶粒和晶粒之间存在着界面，即晶界。晶界处晶粒从一个位相过渡到另一个位相，晶界把结构相同位相不同的亚晶粒隔开。由于晶界处同时受到两侧晶粒的影响，处于两侧晶粒的

25、过度状态，从而使晶界具有特殊的构型，进而对晶体的结合性能产生着巨大的影响。如图1.1所示为Cu晶体中9（221）=141.1°的晶界结构模型13。上下两部分晶粒的晶向不同，交叠在一起，在界面处形成E结构，晶界处为两边晶粒的过度状态。晶界对多晶材料的物理、化学、力学等性质有着非常显著的影响。材料的强度和断裂等力学行为，如偏聚、晶界扩散，以及晶界处初始滑移位错的形核等，都受到晶界结构的显著影响。图1.1 Cu晶体中9（221）=141.1°的晶界结构模型13晶界按晶粒间取向差的大小可分为小角度晶界（取向差小于5度）和大角度晶界（向差大于10度）。小角度晶界通常有可分为扭转晶界和

26、倾侧晶界两种类型，如图1.2和图1.3所示。重合位置点阵模型（CSL）14：即在一些特殊位相的晶界处，有一些原子同属于两边晶粒的格点，且自身形成了超晶格点阵模型。图1.2为面心立方结构中的（001）面重位扭转晶界15，即当旋转角=36.9度时所得的扭转界面，其中重合点阵的格点（即图1.2中的大黑点）的数目相当于总格点数目的1/5，即为重合密度。设为重合密度的倒数。越小，即界面处重合的密度越大，重合的原子数也就越多。图1.2 面心立方结构中的（001）面的扭转晶界15图1.3为 <100>38度时，重位倾侧界面示意图15。图中AC为阶，阶高为BC，阶长为AB。显然，阶越小，即图1.3

27、中AC越短，重合密度就越高，也越小，阶中不接触的原子也越少，即晶界能也较小。图1.3 <100>38°时，倾侧界面示意图15相界为不同两相所形成的界面，相邻两相不仅取向不同，而且结构、成分也不同。按照原子在相界上排列不同，可把相界分为三种形式16：（a）共格相界面：界面两边金属具有同种类型的晶格结构，如都为FCC结构，且其晶格参数相差不大；（b）半共格界面：界面两边金属，其晶格类型相同，但其晶格不匹配程度稍微大点；（c）非共格界面：界面两边金属具有不同的晶格结构，如FCC/BCC界面。通常将失配度定义为一个量，= 2|d1- d2|/(d1 + d2)，这里d1, d2分

28、别晶体1和晶体2中原子间的距离。两近邻平行失配位错间的距离L=b/，这里b为失配位错的Burgers矢量。在共格和半共格晶界或是相界面上，由于晶格不匹配，经弛豫后，在界面处会形成失配位错网络，而这些界面失配位错网络对界面的初始滑移位错的形核和发射起着决定性作用，从而决定着界面的结合强度。金属多层膜之所以具有其都特的性能，主要归因于存在着界面。界面既可以阻碍位错运动，也是位错的发射源17-21。位错和界面相互作用决定着金属多层膜机械性能。关于位错与界面的相互作用机理的研究就变得非常有价值和意义。1.3晶体位错相关理论概述位错是指晶体中滑移部分与未滑移部分的分界线，是晶体中原子的局部不规则排列，为

29、一维缺陷或线缺陷。晶体的塑性变形通过位错滑移来实现，若没有位错的运动，变形就不能发生，位错滑移通常在特定的滑移系上进行。其存在对材料的物理性能、化学性能，以及机械性能产生着巨大的影响。二十世纪初期，Volterra等人22在连续弹性力学中提出位错的概念。实际应用中的材料，其内部通常含有这样或那样的缺陷。缺陷的存在对材料各方面性能会产生巨大的影响。如会降低材料的使用寿命、力学性能等等。但是有时候材料中存在缺陷又会增加其某方面的性能，如在半导体中掺进微量元素来提高其性能，加工增加位错来使金属材料强化等。因此，为了使材料更好的为人类服务，对材料中缺陷的研究就显得非常重要。晶体力学性质的微观理论和位错

30、有着非常密切的联系。晶体的塑性变形通过滑移来完成，在不同的滑移系下进行，而晶体滑移借助滑移位错来实现。近年来，随着晶体力学性质的微观理论和实验的发展，进而促进了晶体缺陷相关理论的发展。根据局部滑移的方式不同，位错可分为刃型位错、螺型位错，以及混合位错，如图1.4所示。图1.4(a)为刃型位错模型，从图中可以看出，有一额外的半原子面终止在晶体中形成一条线缺陷，即位错线，图中用符号表示。由图1.4(a)可以看出位错线上部分原子受到压缩，而下部分受到拉伸，在位错线处形成局部应变场，导致位错线处晶格的对称性遭到破坏，该处相对完整晶体，势能会增加。当额外的半原子面在晶体的下半部分时，用标示。其多出的半原

31、子面如同刀刃一样插入晶体中，故其形成的位错称为刃型位错。 (a) 刃型位错23 (b) 螺型位错24(c) 混合位错24图1.4 三种类型位错示意图图1.4(b)为晶体中螺型位错的晶体外观。晶体上下两部分沿滑移面相互滑移了一个原子距离，滑移终止在晶体中。螺型位错也是一种线缺陷，如图1.4(b)中的AB, 即为位错线，是已滑移区与未滑移区的分界。大部分晶体材料中的位错既不是纯粹的刃型位错也不是纯粹的螺型位错，他们往往是这两种类型都有，即为混合位错，如图1.4(c)所示。通常用柏氏矢量来表示晶体滑移的量。下面以刃型位错为例介绍一下柏氏矢量。在图1.5(a)中封闭回路MNOPQ，而在图1.5(b)中

32、相应的线路是不封闭的，始点M，终点Q。我们定义：将终点联到始点的矢量（图中的），称为柏氏矢量。在图1.4三种类型位错示意图中，其中矢量b即代表柏氏矢量。图1.5 刃型位错中的柏氏矢量示意图151.4晶体界面的分子动力学研究现状在过去十年里，对金属多层膜各方面性能的理解取得了巨大的进展。纳米尺度的金属多层膜在屈服应力、塑性、抗腐蚀性能等方面具有特殊的性能6-11。这种特殊性能主要归因于多层膜中存在着界面，界面既可以阻碍位错运动，也是位错的发射源17-21。获得精确工程界面的能力在当今新技术发展方面扮演着越来越重要的角色25。位错和界面的相互作用对金属多层膜机械性能的影响依赖于界面的类型。过去几十

33、年里，分子动力学方法已被广泛应用于材料科学、物理学、化学、生物学，以及医学等领域。随着高性能计算机的发展，原子模拟在材料性能预测与设计方面已成为一种有效的方法。已经有大量的研究者用分子动力学方法研究晶界和界面处位错的形核和发射过程26-34。如，Tschopp等人35用原子模拟方法研究了纳米晶Cu中位错的形核和发射过程，表明晶界是位错形核和发射源；Spearot等人28 用分子动力学模拟方法研究了单轴拉伸下位错从双晶界面的形核和发射过程；Derlet等人17用原子模拟方法模拟了位错和Cu晶界的相互作用，表明晶界既可以阻碍位错运动，也是新位错的发射源；Shao等人26 用分子动力学模拟研究了Ni

34、/Cu双层膜的纳米压痕过程，结果表明界面会阻碍滑移位错位错穿过界面；Cheng等人36用分子动力学方法研究了Cu/Ni双层膜界面的初始失配位错网络，以及其对摩擦过程的影响，结果表明界面失配位错网络对位错穿过界面有明显的阻碍作用；Zheng等人19用分子动力学方法模拟了Ni/Al双层膜的单轴拉伸过程，表明初始滑移位错从界面的失配位错线形核和发射等。然而，大部分这些研究主要关注于同种类型的材料和非扭转界面，如晶界、FCC/FCC界面。不同类型的材料，如FCC/BCC界面，以及扭转界面被研究的相对较少。对于FCC/BCC界面有如下：Hoagland等人对非共格的Cu/Nb界面进行了研究18,37,其

35、主要关注于在二维条件下位错和界面的相互作用； Wang等人38用类似方法更细致地研究了位错和不同类型界面的相互作用，认为Cu/Nb作为非共格界面系统，界面会阻碍滑移位错的发射；更近一点，Shao等人16用分子动力学方法研究了Cu/Nb双层膜的纳米压痕过程，表明此非共格界面对滑移位错的增殖起到强烈的阻碍作用。而对于扭转界面和扭转晶界的研究也非常少，这其中包括在铜扭转晶界的晶界迁移和晶界自扩散的研究39；未弛豫的Ag(001)/Ni(001)扭转界面能分析40；不同扭转角度下Cu纳米线机械性能的研究41等。这些工作对FCC/BCC界面，以及扭转界面和扭转晶界的研究提供了非常宝贵的信息。本文用分子动

36、力学模拟研究了位错和BCC-Fe/Ni界面的相互作用，以及不同扭转角下Cu(001)/Ni(001)双层膜的结合性能。1.5本文的研究思路、目的及意义薄膜材料在当今社会中具有重要的应用价值，目前它已被广泛应用于航空航天、机械制造、电子技术、光学工程以及计算机工程等各个领域。而在薄膜材料的应用过程中，薄膜的使用寿命和可靠性是人们普遍关注的焦点问题。而界面的结合性能是影响多层膜寿命和可靠性的关键指标。位错和界面的相互作用机理决定着界面的结合性能。即位错和界面的相互作用机理在薄膜寿命和可靠性方面扮演着关键角色。因此对位错和界面的相互作用机理的研究就显得特别有价值和意义。位错和界面的相互作用机理虽然已

37、经有很多科研工作者在研究，但是大部分研究集中在位错和FCC/FCC界面，以及位错和晶界间的相互作用，而位错和FCC/BCC界面，以及位错和扭转界面的相互作用机理研究的相对较少。同时，近年来分子动力学模拟方法已成功应用于位错和界面相互相互作用的研究中。Fe-Ni合金作为一种重要的结构功能材料，由于具有较好的抗腐蚀性能、断裂韧性，以及较好的塑性42,43，而被广泛应用于许多领域。下面简要介绍下本文的研究思路、研究目的和研究意义。研究思路是首先建立Fe/Ni双层膜模型和Cu/Ni扭转界面薄膜模型。然后用分子动力学方法模拟了这两模型的拉伸过程。Fe/Ni双层膜采用侧向拉伸，拉伸应力仅应用于Fe薄膜一侧

38、。Cu/Ni扭转薄膜模型拉伸过程中，沿轴向拉伸。研究了弛豫后界面的初始失配位错网络，以及其对界面位错的形核和发射的影响，进而对界面结合性能的影响，以揭示位错和界面的相互作用机理。研究目的是用计算机模拟方法研究位错从界面失配位错线形核和发射过程的微观机理, 以及如何设计薄膜材料使其界面结合强度最高，使材料更好地为人类服务。本文研究的意义在于用分子模拟来研究金属/金属薄膜材料界面的结合性能，揭示界面位错的微观演化机理，为实验和理论研究金属/金属薄膜材料界面结合性能提供一定的支持和帮助，使之为人们更好的服务。第2章 分子动力学方法2.1 引言如今，随着科学技术的不断发展，人们对自然世界的探索已不局限

39、于宏观领域，越来越趋于微观领域，如分子、原子尺度。由于宏观物质是由无可计数的原子、分子构成，处理起来十分复杂，甚至超越了人的计算能力，以致处理起来十分复杂。直至计算机模拟技术的出现, 是这一问题处理起来方便许多, 另外也降低了成本，从而大大的促进了微观领域的发展。计算模拟是化学、材料学、医学、生物学等多个领域的交叉学科，是一门正在蓬勃发展的新兴学科，是利用一定的算法结合相关计算机模拟软件对物质的组成、微观结构、性能以及服役性能进行计算模拟与设计的学科。分子动力学模拟44是指用计算机模拟原子或分子的运动过程, 以揭示材料内在的本质规律。近年来，分子动力学模拟已成为纳米领域研究最常用的计算方法。与

40、蒙特卡洛方法、第一性原理，以及有限元一起成为计算机模拟技术的主流技术，被称为计算机实验。计算机实验在当今科学研究领域扮演着重要角色。在计算机实验中，通过建立理论模型，并借助计算机以及一定的算法和程序来计算。其中，蒙特卡洛方法是研究随即性问题最有效的方法，可以研究材料介观尺度问题。第一性原理方法可以无需任何实验数据，完全从材料组成的种类以及排列方式出发计算材料性能，其研究体系比分子动力学小得多。有限元则是研究宏观尺寸的有效方法。这些计算机模拟方法被广泛的应用于科学研究，为理论研究和实验研究搭建桥梁，成为除理论研究和实验研究外探索自然界的另一种有效方法。模拟有时是一种理论，但有时又是一种实验。一方

41、面，我们处理的是模型，而不是真实的物体，这表明模拟属于一种理论方法；另一方面，通过计算机模拟证实一个模型的程序又非常类似于实验：我们运行一个程序，然后分析结果，其方法非常类似于实验过程。计算机模拟、理论研究和实验研究三者相互联系，相互补充，相互促进，共同推动着科学技术的发展。计算机模拟可以方便的观察到许多实验无法观察到的细节，甚至实验很难达到或是无法达到的条件，揭示材料的本质规律。鉴于计算机模拟的优点，本文使用分子动力学方法来研究位错和Fe/Ni界面的相互作用机理，以及Cu/Ni界面的结合强度。2.2基本原理分子动力学模拟通过用经典的牛顿运动方程来描述粒子的运动过程。2.2.1 积分方法常用的

42、积分方法有Velet算法、Leap Frog算法、Velocity-Verlet算法和Beeman算法等。(1) Verlet算法在分子动力学模拟中，最常用的时间积分方法为Verlet45 算法。其对粒子的坐标进行泰勒展开得到式2.1和式2.2： (2.1) (2.2)式中r代表位移，v代表速度，a代表加速度，t代表时间。将式2.1和式2.2相加并忽略高次项得到t+t时刻的位置，即式2.3。 (2.3)式2.3即为Verlet算法的基本形式。将式2.1和式2.2相减并忽略高次项可得粒子运动的平均速度 (2.4) 或者 (2.5)Verlet算法具有算法简单，存储要求适度等优点，但容易造成精度损

43、失，且非自启动算法。（2） Leap Frog算法Leap-Frog算法46，即蛙跳算法。其表达式如下： (2.6) 即得到 (2.7) 由式2.5得 (2.8) Leap-frog算法比Verlet 算法计算效率高、计算量小。但其缺点是位置和速度不是同步的。（3）Velocity-Verlet算法Velocity-Verlet算法47其优点是计算量适中，且给出了显式速度项，目前得到了广泛的应用。其表达式如下： (2.9) (2.10) 在Velocity-Verlet算法中前一个时刻的速度、力和位置需要被储存。（4）Beeman算法Beeman算法48也是在Verlet算法基础上发展的另一种

44、计算精度更高算法，但增加了其计算量。2.2.2 原子间的相互作用势在分子动力学模拟中，势函数是影响计算结果的关键因素之一，其决定着计算结果能否准确地反映实际材料的性质。原子间的相互作用力也可由势函数求得，若势函数存在误差，则所计算的相互作用力也会受到影响。势函数的准确性决定着分子动力学模拟的精度。(a) 对势对势通常有Lennard-Jones势、Morse势和Johnson势等势模型。(1) Lennard-Jones（L-J）势49： (2.11)式中为能量参数，rij为原子i和原子j之间的距离，为长度参数。当rij=时，势能为零；当rij时，表现为引力势，由式中第二项来表示，其最小值为1

45、.122, 当距离较远时，原子间的引力会越来越弱，趋近于零，通常为了计算方便，选取截断半径rc, 当rijrc时，取势能为零，在Lennard-Jones势模型中，截断半径rc通常取21/6；当rij时，表现为斥力势，由式中第一项来表示，当原子间的距离越近，部分电子云发生重叠，原子间的斥力会越来越大。(2) Morse势50： (2.12)式中A、r0分别是结合能，调节参数和平衡间距，这三个值都是取经验参数。金属固体常用这种势函数来描述，常用于金属铜。(3) Johnson势51 (2.13)-Fe常用Johnson势函数来描述。(b) 多体相互作用势实际上，在多粒子系统中，就不能简单的应用以

46、上对势模型，而必须考虑多体间的相互作用。常见的多体势有嵌入原子法、Stillinger-Weber势、Finnis-Sinclair势等。(1) 嵌入原子势（EAM）52在嵌入原子势（Embedded Atom Method,简称EMA）中，材料中的每一个原子被看作如同一个杂质原子镶嵌在包含所有其它原子的基体中。EMA势可以处理各种各样的金属体系，其中包括断裂、表面、杂质，以及合金等体系。系统的总能量如下： (2.14)EAM势是将晶体的总势能分成两部分组成，上式中右端第一项ij代表对势相互作用，rij为原子i和原子j之间的距离。第二项Fi即嵌入能，是电子密度i的函数，代表多体相互作用。式中1

47、/2表示对势为两原子共有，式中i可以表达为： (2.15) 对于不同的金属，EAM势需通过拟合金属的宏观参数而得到。(2) Stillinger-Weber势53Stillinger-Weber势常用于描述半导体材料。 (2.16)式中jik 为ij键和ik键之间的夹角，g(r)为截断半径在第一近邻和第二近邻之间的衰变函数。截断半径的选取是为计算方便，因为当距离非常远时，原子间的相互作用势非常小，趋近于零，当截断半径取得越大，结果越精确，但计算量也越大。当cosjik=-1/3时，这种结构类似金刚石四面体结构，其结构非常稳定。(3) Finnis-Sinclair势54以合金属能带紧束缚理论为

48、基础而发展的Finnis-Sinclair势，给出了详细的多体作用势。其表达式如下： (2.17) (2.18)式中，当x>0时，H(x)=1；当x<0时，H(x)=0。Ak、Rk、ak、rk为常数，且有R1>R2，r1>r2>>r6。它们的值随随具体应用的材料而定。2.3 边界条件选取合适的边界条件，对分子动力学模拟非常重要。但通常主要考虑以下两点：首先，为了减小计算量，应尽可能减小模拟的单元，但是为了排除可能的动力学扰动对结果造成影响，以及为了满足统计学处理的可靠性要求，又要求模拟原胞应该足够大；其次，要从物理角度考虑应变相容性、体积变化及环境的应力平衡

49、等实际耦合问题。当模拟大块体时，需使用三维周期性边界条件；当模拟薄膜材料时，需使用二维周期性边界条件，如取X、Y方向为周期性边界条件，Z方向为自由边界条件；当模拟纳米线时，可选取一维周期性边界条件，如取Z轴方向为周期性边界条件，X、Y方向为自由边界条件；有时又需使用混合边界条件，即几种边界条件的组合。实际应用时应根据实际情况选择合适的边界条件。由于实际的宏观物质由无可计数的原子、分子组成，而计算机的计算能力又有限，所以为了进行模拟，就必须采用周期性边界条件。当使用周期性边界条件时，实际的粒子被假想在一个盒子里，这个盒子在空间三维方向被重复无穷多次，填充满整个空间。即如果在盒子中位置r处有一粒子

50、，那在式2.19所在位置处也能找到相同的粒子。 (2.19)这里l, m, n为任意整数，a, b,c 为大小对应盒子边长的矢量。所有这些镜像粒子一起运动，但实际在计算程序中的只有一个粒子。这样大大的减小了计算量，而且也不失精度。图2.1为分子动力学元胞三维周期性边界的二维图示。在基本分子动力学元胞中，当有原子离开分子动力学元胞时就有对应原子进入其中。 图2.1 分子动力学元胞三维周期性边界的二维图示2.4 温度、压力控制方法在分子动力学模拟中，我们往往要考虑模拟系统是否与外界发生作用，如温度、压力等。对于不同的模拟对象，往往要对原子数N、体积V、压强P、温度T等量中的一些进行控制，以实现不同

51、的系综。如在NVT系综下需要控制原子数N、体积V、温度T保持不变。而在NPT系综下则需要对原子数N、压强P、体积V进行控制。2.4.1 控温方法在NVT、NPT系综下，甚至NVE系综下，我们期望模拟体系温度保持恒定不变。而系统的温度T与动能K及速度有直接的关系，如下， (2.20)kB、i分别为Boltzmann常数和第i个原子的速度，N为原子个数， Nc为约束自由度数。温度控制方法通常有速度标度法、Berendsen热浴法和Nosé-Hoover控温方法等。(1) 速度标度法55由式2.20可知温度和速度存在着一定的对应关系，故可通过对速度进行标度进而达到控制温度的目的。若速度乘以

52、标度因子后，t时刻的温度的变化为 (2.21) (2.22)Treq，分别为期望的参考温度和速度标度因子，给每一步乘上标度因子来达到控制温度的目的。(2) Berendsen热浴法56Berendsen热浴法是假想系统和一热浴相接触。对速度每一步都进行标度，标度因子为 (2.23)式中，t，tT 分别为时间步长和耦合参数，Tbath，T(t)分别为t时刻系统的期望温度和瞬时温度。（3）Nosé-Hoover控温方法57-58Nosé-Hoover控温方法中，一恒定温度的热浴与体系相耦合。2.4.2 控压方法在NPT, NPH系综下，我们需要模拟定压情况下的分子动力学问题，下

53、面主要介绍Berendsen压浴法和Andersen方法。（1）Berendsen压浴法59Berendsen压浴法是通过对体积乘以标度因子Cp，来通过控制系统体积来实现对压力的控制。如同让一压力为期望值的压浴与系统接触。 (2.24) 式中，P(t) 和Pbath分别为系统在t时刻的瞬时压力和期望压力；k和tP为耦合参数。(2) Andersen方法60Andersen方法类似于图2.2中的活塞调节压力，通过活塞调节体积变化来控制压力。图2.2 活塞法调节压力2.5 分析方法2.5.1 中心对称参数（centrosymmetry parameter）为了可视化晶体缺陷，特别是位错的形核和发射过程，我们需要在晶体中可靠的确定这种缺陷结构的位置。基于不同的标准有许多不同的方法来研究晶体缺陷，如原子的势能、配位数，以及中心对称参数等。本文采用Kelchner等人61介绍的中心对称参数方法。在固态物质中，中心对称参数对于区分原子为完整晶体，局部缺陷（位错或堆垛层错）或表面原子非常有用。其定义为： (2.25)这里 为近邻原子对数。对于 BCC结构，等于4；FCC结构，为6。Ri 和Ri