高介电电子材料Hf02合金的热学性质模拟

1、毕业设计（论文）高介电电子材料Hf02合金的热学性质模拟 学 生 姓 名 xxx 学 号 07080401 专 业 班 级 高新材料基地二班指 导 教 师 xxxx 提 交 日 期 xxxxx 材料科学与工程学院51摘 要主本文首先通过计算机模拟所得到HfO2氧化物合金材料的局域晶格结构和HfO2氧化物合金材料的原子振动情况。根据原子间相互作用力等微观信息了解多原子体系的结构和性质等，采用分子动力学方法，就是数值地分析求解多体系统的确定性运动方程，并根据对所求结果进行统计处理决定粒子的轨迹，从而统计出HfO2在不同温度下融化表面原子发生的运动变化，初步画出加入不同温度，加热不同时间后，晶格八个

2、层面上的原子的均方位移情况，并绘制出关系图，再通过比较得出系统的曲线对比结论。关键词：扩散，振动，表面结构，分子动力学模拟Abstract The paper first received by computer simulation of HfO2 - SiO2 oxide alloy material of local grid structure and HfO2 - SiO2 oxide alloy material of atomic vibration condition. According to the atomic interaction force more micro i

3、nformation about such structures and properties of the atomic system, using the molecular dynamics method, the numerical solution is to analyze the uncertainty of multi-body system, and according to the motion equation of statistical processing results of particle trajectories, Thus statistics a HfO

4、2 - SiO2 at various temperatures melt surface atomic happen movement change, preliminary draw to join different temperature, heating different time, the lattice eight level of atoms, and both azimuth shift mapped relation graph, again through comparison of the curve that system conclusions. KEY WORD

5、S: Diffusion，vibration, Surface structure, Molecular dynamics simulation目 录摘 要1Abstract2目 录3第1章 绪论41.1 课题研究背景41.2 课题的研究目的及研究内容51.2.1 研究目的51.2.2 研究内容5第2章 原子模拟方法简介62.1 原子模拟方法62.2. 晶体模拟中的基本原理72.3 计算晶格能的方法102.4 动态模拟142.5 计算程序152.6 原子模拟成果举例15第3章 分子动力学简介与GULP软件163.1 分子动力学简介163.2 分子动力学（Molecular Dynamics,简

6、称MD）方法163.2.1 引言163.3 分子动力学基础知识183.3.1 基本理论183.3.2 原子间相互作用势213.4 分子动力学在材料科学中的应用223.4.1分子动力学的概述223.4.2分子动力学的适用范围223.4.3分子动力学的应用233.4.4分子动力学的最新进展253.5 Gulp软件在分子动力学模拟中的应用253.6 运用GULP软件进行均方位移(MSD)数据分析27第4章 模拟结果分析与讨论284.1.均方位移(MSD)分析324.1.1均方位移MSD324.1.2 A(110)面原子在不同温度下运动方式随时间推移的变化曲线33参考文献37外文翻译39中文翻译48致

7、 谢54第1章 绪论1.1 课题研究背景随着半导体集成电路制造工艺的不断升级1,2，当传统栅介质材料SiO2的厚度减小到纳米尺寸时，电子将因量子隧穿效应直接穿过介质层，导致器件失效。为了解决这一问题，采用高介电常数（High-k）的材料来替代SiO2形成栅介质层是一个重要的技术途径。若欲取代SiO2成为半导体器件中的栅介质，High-k材料必须具有与SiO2/Si系统相似的晶格和电子性质，并要与当前的半导体制造工艺相兼容。介电值材料常以High-k称呼之，主要功能是用来隔绝闸极的漏电流，High-k意指高介电常数，是用以衡量一种材料能储存多少电荷。High-k材料，当k值越高，晶体管的电容值也

8、越高，也就能降低晶体管温度，控制漏电流，且High-k介电值厚度是二氧化矽之数倍，能够有效降低漏电流达100倍3Yw9D S m6d8VKa#。High-k制程在半导体业界并不是新发现，只是制程的应用在2007年才逐渐成熟，Intel和IBM在2007年1月相继发表45nm时代的实用化时程，也就是所谓的High-k/metal gate材料技术。Intel已在Penryn的芯片上，应用High k/ Metal Gat未来将扩往32nm线宽的晶体管IC。而IBM也成功应用High k/ Metal Gat在32nm的SRAM芯片。总体而言，High-k材料的实验研究侧重于继续寻找新的氧化物或氧

9、化物“合金”，寻求更好的电学性能和制备/生长工艺；理论研究多侧重于利用第一性原理探讨原子分布和电子性质。但实验研究往往难以得到原子层次上的细节，而第一性原理计算量太大只能处理较小的体系（不多于200个原子）。由于目前集成电路仍以Si为衬底，考虑到High-k材料与Si的晶格匹配和工艺兼容，我们拟采用原子模拟技术来考察以SiO2为基础的High-k材料（HfO2/SiO2），调研发现该材料在原子层次上的信息还比较缺少。原子模拟方法计算速度快，我们将在晶格尺寸原子层次上在较大体系（几千到几万个原子）的基础上详细考察该材料的原子分布、局域晶格结构和它们与介电性能的相关关系。1.2 课题的研究目的及研

10、究内容1.2.1 研究目的1. 给出HfO2氧化物合金材料的局域晶格结构；采用分子动力学的模拟方法研究加热下HfO2纳米团簇的熔化与扩散行为，主要采用的统计方法有Z方向的概率分布函数、径向分布函数、结构有序参数、X射线衍射强度、均方差或均方位移。采用C语言编写了相关的程序用于统计HfO2纳米团簇的模拟数据。2. 得到HfO2氧化物合金材料的原子振动情况。着重研究原子由于热运动将偏离其理想晶格位置，该位移的均方根随时间的变化可以反映出原子在做振动还是扩散运动。原子团簇在加热初期，其均方根位移几乎不随时间增长而增大，说明原子处于振动态；而到了后期其均方根位移随时间增长持续增大，说明该原子进行了扩散

11、运动。1.2.2 研究内容1. HfO2氧化物合金材料的局域晶格结构随温度压力等变化，随着熔化时间的增加，原子团簇中原子的位移距离发生了很明显的变化，在加热条件下原子运动的无序化程度随着时间增加而不断地变大。2. HfO2氧化物合金材料的原子振动随温度的变化，主要根据原子间相互作用力等微观信息了解多原子体系的结构和性质等，所采用的都是基于统计理论的数学解析法。采用分子动力学方法，就是数值地分析求解多体系统的确定性运动方程，并根据对所求结果进行统计处理决定粒子的轨迹，从而给出物性预测和微观结构信息的一种模拟方法，统计出表面原子在不同温度下运动轨迹随时间变化而发生的振动或者扩散现象。第2章 原子模

12、拟方法简介2.1 原子模拟方法经典原子模拟最早是由Boswara和Lidiard进行的，他们试图确定NaCl结构和卤化铈中肖脱基缺陷（Schottky defect）的形成能8,9。实际上，早期的模拟绝大多数考虑高离子化和极其简单的化合物。在20世纪70年代用类似的方法研究了过渡金属氧化物10,11。Catlow研究组的工作是原子模拟技术发展的主要推动力，早期集中于UO2缺陷能以及核裂变产物的计算。在这些研究中，很多先驱性的技术得以应用。原子模拟既可以进行完美晶体的计算，也可以计算缺陷能。开始时原子模拟仅限于块体状缺陷性能的计算，现在已经能够计算表面性能了。计算机工业在过去的几十年里快速发展，

13、集成电路上晶体管的数目一直遵循摩尔定律按指数方式增长，计算机计算的功能越来越强大。最近，计算机原子模拟已应用于生物领域，用来确定复杂的蛋白质结构。这在十年前是不可想象的复杂系统。原子间的力可以通过两种方式来确定：基于牛顿力学的经典方式和基于量子力学的第一性原理的方式。虽然我们采用经典方法用计算机来模拟，但不应该忽视量子力学方法的研究，因为从本质上看后者更为精确。但量子力学方法计算量巨大，当面对有限的计算资源时，它就不太适应性了。对于原子模拟的每个离子，经典计算远为高效，且计算量要小的多，这使得它能够考虑大量数目的原子。当然，当未来计算能力有了根本性的飞跃之后，对更大的系统用量子力学模拟将是更好

14、的选择。有些研究者对晶格的局部结构采用量子力学方法，对整个晶体采用原子模拟技术，不失为一种较好的折中方法。在这一章，我们简要介绍原子模拟技术的基本原理和方法：离子晶体中，单个离子用核壳（core-shell）模型描述12，离子间有长程（库仑作用）和短程相互作用（本论文采用Buckgham形式），采用爱瓦尔德（Ewald）方法计算晶格能，然后使用两种优化技术得到0K和一定温度下的最小晶格能，晶格结构随之确定。还说明了该技术适用的范围和取得的成就。最后介绍我们采用的掺杂LaMnO3和La2CuO4的短程势参数。还考察LaMnO3势参数在压力下和升温后（Ca/Sr掺杂）的情况以进一步检验这些势参数。

15、2.2. 晶体模拟中的基本原理铁电材料在晶格场的作用下，离子的偶极矩的大小和方向很难保持不变，因此离子晶体的模拟常需要考虑极化现象，尤其是当离子半径比较大和有晶格缺陷存在时。在原子模拟的模型中，极化以两种方式存在。在第一种方式中，离子因为承受临近离子的作用而稍稍偏离它们的平衡位置，这称为位移极化。在第二种方式中，是绕固定离子的电子（密度）对原子核的位移，这称为电子极化。在模拟方法中加进电子极化是很重要的，尤其考虑大半径的离子和有电荷缺陷时。然而，包含电子极化增加了自由度，因此也增加了计算量。为考虑电子极化，原子模拟技术采用了Dick和Overhuaser12所设计的核壳模型（图2.1）。这个模

16、型认为每个离子是由一个带X电荷（单位为e）的核和一个带Y电荷的对应壳构成，因此，离子的总电荷是X+Y。无质量的壳与包含离子全部质量的核通过一个弹性力常数k耦合，由此，自由离子的极化率，可以表达为：（2.1）其中是0真空的介电常数，如果Y的电荷单位采用电子电荷，则k的单位为eV·Å2。参数Y和k可以用来拟合介电常数和弹性常数。这个拟合主要考虑了高频介电常数，因为它来自于电子极化，而不是像静态介电常数那样来自离子极化的贡献。尽管核壳模型在本质上是唯象的，但它取得了广泛的成功。与其它模型如硬球离子模型（rigid ion model）或点极化离子模型（point polariza

17、ble ion model）相比，这个模型的优点在于，作用于一个离子上的任何力都被认为是通过壳来作用的，因此与短程交互作用相联合而导致极化，见图2.1。显然，核壳模型提供了一个框架，与其它模型（如硬球模型）相比，核壳模型可以模拟离子间更多的交互作用。XY 图2-1 核壳模型：实心小圆代表离子的核，虚的大圆和它上面的空心小圆代表壳，双向箭头表示离子自身的极化，实直线表示库仑作用，虚线表示短程作用。离子晶体（例如氧化物）中离子间存在着两类相互作用。一类是长程相互作用，就是异种电荷离子间的库仑引力和同种电荷离子间的库仑斥力。这类离子间的作用力有特别简单的形式，就是库仑定律：（2.2）其中和分别是离子

18、i和j的电荷，是离子间距，是真空的介电常数。另一类是短程相互作用。两个相邻原子如果靠得足够近时，其电荷分布能够重叠。这导致排斥的相互作用，如果这些原子间的距离变的足够的小，它们之间的总力则为斥力，即使这些离子带异种电荷。排斥是由两个方面带来的，其一是泡利项，这是泡利不相容原理的结果。其二是原子核之间的相互排斥。短程排斥力的作用是：虽然异类带电离子相互吸引，但短程作用防止它们吸引得太近而崩塌；虽然同类带电离子相互排斥，但短程作用防止它们排斥得太远而解离。当核间距较小时（但比前面斥力所需的间距大），还存在一个引力，范德华伦敦（Waals-London）交互作用。这是一个相对弱的力，它来自于每个相互

19、作用离子之间产生的自发诱导极化19，像德拜（Debye）猜想的那样。利用量子力学的方法，伦敦（London）确定了这个引力的一个通用表达式。由于电子的关联运动产生了偶极子20，在两个相同原子的情况下，这个力与呈线性关系。尽管这是一个量子力学效应，线性依赖关系仍可以从经典静电学中得出。对两个H原子之间范德华伦敦力的完整推导可在文献21,22中找到。然而，应该注意的是，这个力与离子的极化有关。基于此，通常不考虑阳离子阳离子的范德华伦敦交互作用，因为阳离子半径趋向更小并且不极化。库仑项与短程项的结合，是由玻恩（Born）和Landé最早引入23的，它可以表示为：（2.3）其中b和n是常数，

20、用以表达平衡态离子间的距离，r是异类离子间的最近距离。采用这个模型的早期工作都取n=9。尽管这是一个极好的近似，尤其是对离子化合物，如碱土金属的卤化物，但通过量子力学计算发现，它并不严格正确，因此后来对这个模型进行了扩展修正。玻恩和Mayer引入了一个短程排斥势24，形式如下：（2.4）其中A和都是可调节的参数。式2.7与式2.8的不同在于式2.8包含了一个指数排斥项。在这一点上，玻恩和Mayer也添加了一个引力项来说明范德华相互作用。与范德华、伦敦和Margenau的工作相一致，这一项也采用了的形式，C是一个可调节的参数。我们的原子模拟研究只考虑固体，因此需要对Lennard-Jones势作

21、一个替换，如果将式2.4的短程斥力项与范德华引力项相结合，可以得出式2.5所示的Buckgham势。（2.5）在Buckgham势中，A、和C都是可调节的参数，它们对短程交互作用的描述很大程度上决定模拟是否能取得成功。图2.2描绘了长程势和短程势对整体离子交互作用的影响，用核壳模型描述的两个离子间的长程与短程相互作用可参见图2-2。如何确定这些参数将在下节里讨论。长程作用总的作用短程作用r (Å)E (eV)0 图2-2表示的短程势的参数可被认为具有较弱的物理意义。例如A与离子的“硬度”有关，描述了离子的尺寸，而C被用来描述范德华交互作用。对交互作用相同的离子，C可以表示为：（2.6

22、）对交互作用不同的离子，C可以表示为：（2.7）其中表示静态偶极子极化。K为有效电子数（即导致极化的电子数）。2.3 计算晶格能的方法给定了离子之间的相互作用，可求得晶体的结合能（晶格能）。本文中所考察的所有化合物都是离子晶体，即它们都是正电荷金属离子和负电荷氧离子的规则排列。并假定所有的离子都是球形带电的，那么这些材料的晶格可以用经典方式来描述25。从原则上讲，固体的内能应考虑所有电子和原子核的位置和动量，但在经典方式中，晶格能可以表示为：（2.8）其中第一项是库仑作用（是真空的介电常数），第二项是短程作用（我们采用式2.10的形式）。一般地，式2.8只考虑每两个离子之间的相互作用，而忽略三

23、个离子及其更多数目离子之间的相互作用，对于大多数系统而言，这种近似已足够了。式2.8中的第一项是库仑能，是离子间的远程相互作用，也是主要的相互作用。异性带电离子交替排列，互相吸引，组成了离子晶体结合能的主要部分，它往往占到90%以上的比例。式2.8中的第二项，是总的短程作用能。温度为0K时的计算具有“静态”特征，即它们没有具体考虑晶格振动或位型熵。在计算式2.8中所给出的晶格能时，库仑相互作用采用数学方法求和，而短程作用，在几个晶格间距内被求和，几个晶格间距之外就被忽略不计。尽管在晶格能公式2.8中库仑项形式上很简单，但实际上它却难以计算，因为库仑力是一个长程力。长程力的空间相互作用的衰减比慢

24、（d为系统的维度）。原子模拟时，长程交互作用带来了的问题相当严重，因为它们能跨越模拟晶胞的一半距离。因此，必须引进一个解决这个问题的方法。使用最广泛的方案，包括我们使用的原子模拟技术，都采用了爱瓦尔德方法来求和。爱瓦尔德设计了一个方法，对一个离子和它所在的周期性晶体场之间的相互作用求和。下面给出这个处理方法的简化形式，仅侧重于描述方法的几个关键方面。这里所采用的晶格由球形互不重叠的离子组成，这些离子具有相等电量的正电荷或负电荷。在一个特定的晶格阵点上总的势可以分为两个分势，表达为：（2.9）其中是倒空间的势，是正空间的势。势是具有点电荷的晶格，其等量但符号相反的电荷以高斯分布（Gaussian

25、 distribution）的方式分布于晶格上。势是具有高斯分布电荷的晶格，固定在每个阵点，与晶格电荷的符号相同。当势加到，总的势减小到点电荷原先的值。把总的势分成两部分的意义在于，可以用一个能够被优化的参数来确定高斯峰的宽度，这样，这两个部分快速并独立收敛。Catlow和Norgett为这个宽度确定了一个最佳值25：（2.10）其中，为高斯峰的宽度参数，N是离子种类的数，V为单胞体积。电荷的分布由和两个部分组成，总的势与宽度系数无关。然而，收敛的速度则依赖于这个参数。马德龙常数的定义说明参考点的电荷分布对它本身的势没有贡献，换句话说，离子感觉不到它们自身的静电势。因此，可被表示为差值：（2.

26、11）其中，是晶格上与实际晶格符号相同的一系列高斯分布势，是这个离子上与实际晶格符号相反的一个高斯分布的势。和与其相关的电荷密度，可以通过傅立叶级数分别表达为：（2.12）和：（2.13）其中，和为系数，是倒易晶格矢量的倍。当和减小，而增加时级数收敛。电荷密度可以通过泊松方程（Possions equation）与静电势联系起来：（2.14）这个关系可以用来确定的表达式：（2.15）其中为单胞体积。当=0时，势趋于无穷。然而，由于假定中性单胞的总电荷为0，因此=0时可忽略。组成项，是以高斯几率分布在参考点的势：（2.16）由式2.16得到为：（2.17）总势的另一项为，在参考点估算：（2.18

27、）在每个参考点，它由三部分贡献组成：与离子j有关的点电荷，在格点j处半径以内包含的高斯分布方式的电荷，以及半径外高斯分布的电荷。确定和后，我们就可以由式2.14得到总的势 （2.19）我们的模拟技术是基于从Born模型归纳出来的壳模型的广泛成功使用，运用这种模型，晶格能E可以被述为： （2.20） 这里的第一个术语是Coulombic能量，被理解为长程有效的电荷间的相互作用，第二个术语是短程作用，其作用被Buckingham势能来表述：, （2.21） 这里的, 和 是拟合参数，为了描述单个离子，与其所依赖的原子区域环境的极化，通常用核壳模型。对于任何离子的核与壳之间，弹性常数K之间的关系用下

28、式可以表示： （2.22） 这里的di是离子i的核与壳之间的相对位移。单个带电荷Y的壳与带电荷为X的核（X+Y是离子的化合价）的极化可以下面这个公式表示： （2.23） 这里Y是壳的化合价，与介电常数有关，K是壳与核之间的力学常数，与声频有关，不管是Y还是K都是拟合常数。 能量最小化方法为使前面的模型能够用于预测完整晶格性能，必须结合晶格能最小化方法使系统达到力学平衡状态26。在本论文中，计算了所有的离子交互作用，每个离子移动一个距离，这个距离在总力场的方向上正比于作用于离子上的力。有两个程序来进行晶格能最小化：恒压过程或恒体积过程。恒体积最小化过程通过离子的坐标决定最小能量，只考虑单个离子的

29、应变。对恒压技术，不只是要通过调整离子坐标，也要考虑到单个离子以及各个方向上的应变来调整晶胞尺度，通过这两个手段来确定最小能量。本论文的模拟结果都是在恒压条件下得到的。2.4 动态模拟上述的模拟方法都是在温度0K条件下进行的，为了检验模拟所采用核壳模型短程势是否稳定或者考察稍高温度下的晶体性质，常需要模拟一定温度下的晶格。这种方法假定晶格的振动频率各自独立，各个振动是随体积变化的量子化振子。为了使晶格能最小，必须消除原子扰动产生的动态压强。动态压强是自由能F对体积V的导数，而特定体积下的自由能可以直接从声子频率获得。材料的晶格模型在一定的温度和压强下可以减小其自由能。我们的途径是调节晶胞的体积

30、和原子位置直到压力变为0，这时压强P是由自由能F与晶胞体积V来决定，适合于立方结构材料：. （2.24）在一定体积计算自由能，然后对其作一个小小调整再重新计算，这时dV就决定着压强。2.5 计算程序本论文中原子模拟用到的计算程序有GULP和PARAPOCS。所用到的势参数是用GULP拟合的，并用GULP产生PARAPOCS的输入文件。利用PARAPOCS可以将晶体学单胞沿三个晶轴长大成为包含合适数目原子的超晶胞。如研究Ca掺杂LaMnO3中的电荷有序时，采用了3×1×3的超晶胞（也就是说，晶胞沿a和c方向长大了2倍，沿b方向不长大）。在超晶胞的基础上，就可以用PARAPOC

31、S程序进行加压、升温和掺杂模拟了。本论文的工作都是用PARAPOCS计算完成的。另外，CASCADE程序也是原子模拟技术的重要组成部分（虽然本论文工作没有用到），它主要用来计算晶格的缺陷能。2.6 原子模拟成果举例原子模拟从上个世纪70年代开始应用，到现在已经取得了很丰富的研究成果。利用这些技术，人们研究了一些离子化合物，主要还是氧化物，如KTaO3、Al2O3、TiO2 、石榴石（garnet）、BaO、MgO和LiNbO3中的缺陷。研究了MgO和CaCO3表面上水的行为，NiO/MgO的晶界和方石英（cristobalite）的可能结构。考察了SiO2/MgO/TiO2的振动和弹性结构性能

32、和TiO2 纳米颗粒的结构稳定性。过渡金属氧化物更是研究的一个重点。研究了LaMnO3和Ba4NaRu3O12的结构。LaMnO3掺杂后的表面、缺陷和压力效应也被研究。研究了LaMnO3和LaCoO3的结构和氧空位的扩散。我们研究组研究了超导体YBaCuO、La2CuO4、Ba1-xKxBiO3、HgBaCaCuO和MgB2 ，以及氧化物LaMnO3和SiO2。第3章 分子动力学简介与GULP软件3.1 分子动力学简介分子动力学方法是一种计算机模拟实验方法，是研究凝聚态系统的有力工具。该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中各种微观细节。它是对理论计算和实验的有力补充。广泛

33、应用于材料科学、生物物理和药物设计等。经典MD模拟，其系统规模在一般的计算机上也可达到数万个原子，模拟时间为纳秒量级。2006年进行了三千二百亿个原子的模拟（IBM lueGene/L）。    分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。在忽略核子的量子效应和Born-Oppenheimer绝热近似下，分子动力学的这一种假设是可行的。所谓绝热近似也就是要求在分子动力学过程中的每一瞬间电子都处于原子结构的基态。要进行分子动力学模拟就必须知道原子间的相互作用势。

34、 在分子动力学模拟中，我们一般采用经验势来代替原子间的相互作用势，如Lennard-Jones势 、Mores势、EAM原子嵌入势、F-S多体势。然而采用经验势必然丢失了局域电子结构之间存在的强相关作用信息，即不能得到原子动力学过程中的电子性质。3.2 分子动力学（Molecular Dynamics,简称MD）方法 3.2.1 引言对于一个多粒子体系的实验观测物理量的数值可以由总的平均得到。但是由于实验体系又非常大，我们不可能计算求得所有涉及到的态的物理量数值的总平均。按照产生位形变化的方法，我们有两类方法对有限的一系列态的物理量做统计平均： 第一类是随机模拟方法。它是实现Gibbs的统计力

35、学途径。在此方法中，体系位形的转变是通过马尔科夫(Markov)过程，由随机性的演化引起的。这里的马尔科夫过程相当于是内禀动力学在概率方面的对应物。该方法可以被用到没有任何内禀动力学模型体系的模拟上。随机模拟方法计算的程序简单，占内存少，但是该方法难于处理非平衡态的问题。另一类为确定性模拟方法，即统计物理中的所谓分子动力学方法（Molecular Dynamics Method）。这种方法广泛地用于研究经典的多粒子体系的研究中。该方法是按该体系内部的内禀动力学规律来计算并确定位形的转变。它首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的

36、坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性, 从而得到系统的宏观性质。因此，分子动力学模拟方法可以看作是体系在一段时间内的发展过程的模拟。在这样的处理过程中我们可以看出：分子动力学方法中不存在任何随机因素。系统的动力学机制决定运动方程的形式： 在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序

37、较复杂，计算量大，占内存也大。适用范围广泛： 原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。 自五十年代中期开始，分子动力学方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算机模拟的重要方法。应用分子动力学方法取得了许多重要成果，例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程的研究。其应用已从化学反应、生物学的蛋白质到重离子碰撞等广泛的学科研究领域。实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制

38、：一个是有限观测时间的限制；另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数目趋于无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。数值求解时的离散化方法： 对体系的分子运动方程组采用计算机进行数值求解时，需要将运动方程离散化为有限差分方程。常用的求解方法有欧拉法、龙格-库塔法、辛普生法等。数值计算的误差阶数显然取决于所采用的数值求解方法的近似阶数。原则上，只要计算机计算速度足够大，内存足够多，我们可以使计算误差足

39、够小。3.3 分子动力学基础知识 3.3.1 基本理论分子动力学方法的出发点是对物理系统的确定的微观描述。这个系统可以是一个少体系统，也可以是一个多体系统。这种描述可以是哈密顿描述或拉格朗日描述，也可以是直接用牛顿运动方程表示的描述。在前两种情况下，运动方程必须应用熟知的表述形式导出。顾名思义，分子动力学方法是用运动方程来计算系统的性质，结果得到的既有系统的静态特性，也有动态特性。而蒙特卡罗方法只能得到系统的位形特性，虽然它也有一个动力学解释。 在计算机模拟出现以前，人们根据原子间相互作用力等微观信息了解多原子体系的结构和性质等，所采用的是基于统计理论的数学解析法。然而，原子相互作用力稍微复杂

40、一些，不用说求解统计力学方程严格解，就是进行数值求解也是一件很困难的事情。分子动力学方法就是数值地分析求解多体系统的确定性运动方程，并根据对所求结果进行统计处理决定粒子的轨迹，从而给出物性预测和微观结构信息的一种模拟方法。一个完整的分子动力学模拟程序主要由以下几部分组成10：(1)读入指定运算条件的参数，如初始温度、粒子数、密度和时间等。(2)体系初始化，即选定初始位置和速度。(3)计算作用于所有粒子上的力。这是几乎所有的分子动力学模拟中最费时的部分。如果不使用技巧的话，力的计算时间与N2成正比。采用几种有效的方法可以用来加速短程力和长程力的计算，这时计算时间与N成正比。(4)求解牛顿运动方程

41、。这一步和上一步构成了模拟的核心。重复这两步直至计算体系的演化到指定的时间长度。(5)中心循环完成之后，计算并打印测定量的平均值，模拟结束。对于不同的系综，粒子所遵循的运动规律不尽相同，因此，使用什么样的分子运动方程需要具体问题具体分析。对于一个复杂的系综，单个粒子所遵循的运动规律已经完全确定，多种类型的粒子间相互作用的力的形式也完全确定，这样对于一个粒子系综内部粒子的运动的描述具有明确的表达公式，但是其运动方程的解析却是不可能求得的。对于一小块样品，其内部粒子数也是很多的，要想精确求出系综中粒子间相互作用和运动状态是不可能的，只有通过少量粒子的运动性质计算来求出宏观性质。由于粒子数过少，在构

42、建的计算胞边缘的粒子较多，其受力状态与内部粒子不同，为了解决这一问题，分子动力学方法需要对计算胞设置适当的边界条件，主要有四种11：(1)自由表面边界条件。计算胞中一般只包含几千个原子，对应10-20摩尔量级的材料。因此，自由表面边界条件适用于研究大的自由分子。(2)刚性边界条件。计算胞的外围是一个具有晶体结构的原子位置被固定的“包层”。“包层”的厚度应大于原子间的作用距离，代表着与计算胞内原子有相互作用的宏观晶体。(3)柔性边界条件。计算胞外围的原子在力作用下，可以有小的位移。研究扩展缺陷如位错运动时，常采用此种边界条件。(4)周期性边界条件。一般用于模拟较大的系统，即在计算胞中，最左（或前

43、、上）面的原子与最右（或后、下）面的原子之间有相互作用。分子动力学方法的基本方程都是线性或非线性的二阶常微分方程，对此人们很早就提出了阿达姆斯(Adams)方法、欧拉(Euler)折线法、改进的欧拉法及其龙格-库塔(Runge-Kutta)方法等，并在实际中得到了广泛应用。就分子动力学而言，主要使用Verlet法和Gear法进行数值积分求解，具体求解方法可参考文献12。所谓分子动力学模拟，是指对于原子核和电子所构成的多体系统，用计算机模拟原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在全部其它原子核核电子所提供的经验势场作用下按牛顿定律运动。图3-1 分子动力学的基本原理流

44、程分子动力学模拟的基本原理：分子动力学将连续介质看成由N个原子或分子组成的粒子系统,各粒子之间的作用力可以通过量子力学势能函数求导得出,忽略量子效应后,运用经典牛顿力学建立系统粒子运动数学模型,通过数值求解得到粒子在相空间的运动轨迹,然后由统计物理学原理得出该系统相应的宏观动态、静态特性。图1所示是MD 模拟过程。MD具体的做法是计算机上求运动方程的数值解。通过适当的格式对方程进行近似,使之适于在计算机求数值解。从使用连续变数和微分算符的描述过渡到使用离散变数和有限差分算符的描述,显然会有误差,误差的阶数取决于具体的近似机制,即所用的算法。模拟首先是规定初始条件。为了确定起见,可令初始位置在格

45、子的格点上,而初始速度则从波尔兹曼分布得出。一个按上述办法建立的系统不会具有所要的能量,而且,很可能这个状态并不对应于一个平衡态。为了推动系统到达平衡,需要一个趋衡阶段。可以通过增加或从系统中移走能量,对运动方程向前积分若干时间步,使系统弛豫到平衡态。接着是物理量的计算阶段,沿着系统在相空间中的轨道计算一切令人感兴趣的量。 模拟中,MD采用周期边界条件和最小镜像原理,可以大幅度减少计算工作量14。周期边界条件是将一定数量的粒子N集中在一定的容积V中,这个容积V称为原胞,原胞周围的部分可以看作是原胞的复制,它们称作镜像细胞。这些镜像细胞的尺寸和形状与原胞完全相同,并且每个镜像细胞所包含的N个粒子

46、是原胞中粒子的镜像,原胞在各个方向周期复制便形成了宏观物质样本。这样只需根据原胞周围的边界条件计算原胞内粒子的运动,幅度减少了工作量。原子间作用势能模型的构造对于MD法的应用至关重要。最简单的偶势模型只考虑两体作用,而与其它原子无关,在模拟中运算量小。20世纪80年代以来,各种经验或半经验的多体势模型迅速发展,特别是镶嵌原子(EAM)17既克服了偶势的缺陷,又不会使计算量太大。3.3.2 原子间相互作用势 在分子动力学模拟中，能够正确的描述原子之间的相互作用势是至关重要的，它决定了模拟结果的可靠性13，因此，分子动力学模拟的很大部分工作被放在了对原子间相互作用势的拟合上。将原子间的相互作用用一

47、个简单的势函数表示，由势函数出发即可以得到物质的各种静态物理性质，如结合能、点阵常数和弹性模量等，这样就把问题大大简化了。通过原子间相互作用势对位移求导数就可以用来描述原子之间的相互作用力，根据前面提到的求解诸多粒子经典动力学方程完成分子动力学模拟，就可以获得物质的一些动态性质，如相稳定性等。基于量子力学的第一性原理计算原子间的相互作用势计算量大，对计算机硬件要求过高，因此在计算机模拟中，人们发展出了各种经验的和半经验的作用势进行使用。最先使用的是对势，后来在对势的基础上又发展出了几种多体势。势参数的拟合在选取适当的势函数形式，并指定F(),f(r)和(r)的函数形式后，我们需要针对自己所研究

48、的金属系统拟合势参数，得到合适的多体势。首先给定势参数的初始值，其次利用该多体势导出该系统的物理性质，例如晶格常数，弹性模量等，并与实验值相比较，然后微调势参数，直到获得满意的结果。通常，还通过拟合的多体势计算能量曲线，与Rose方程进行比较以获得对势参数的进一步认可。对于一个由A、B 组成二元合金系统，一个完整的多体势应有三组势参数：即描述同种原子A-A、B-B 之间的作用势的势参数以及描述异种原子A-B之间的作用势（也称为交叉势）的势参数。因此，需要具有纯金属A、B以及AB所形成的化合物的一些物理性质数据用于势参数的确定。对于纯金属以及生成热为负的合金系统的化合物其数据一般可以通过工具书查

49、找，但对于一些生成热为正的系统，由于化合物不能稳定存在，其物理性质数据难以得到。为了解决这个问题，研究组尝试利用第一性原理辅助构建多体势来解决这个问题，通过第一性原理计算得到合适的亚稳化合物及其物理性质数据，用于交叉势势参数的拟合，目前已经取得较好的效果。前面获得的多体势参数是通过对物质静态性能拟合得到的，为了确保多体势的准确性，还需通过分子动力学模拟验证一些动态性能，如晶胞稳定性、相能量序、熔点等。如果利用构建的多体势模拟出来的动态性能与实际出入太大，则需要继续寻找新的势参数，直到静态性能和动态性能都与实际吻合较好，因此多体势的构建是一个极度耗时的实验校正（trial and error）过

50、程。获得了一组优良的多体势参数后，就可以利用此多体势进行分子动力学模拟，得到自己所需要的结果。3.4 分子动力学在材料科学中的应用3.4.1分子动力学的概述分子动力学(Molecular Dynamics,简MD) 用于计算以固体、液体、气体为模型的单个分子运动,它是探索各种现象本质和某些新规律的一种强有力的计算机模拟方法,具有沟通宏观特性与微观结构的作用,对于许多在理论分析和实验观察上难以理解的现象可以做出一定的解释14。MD方法不要求模型过分简化,可以基于分子(原子、离子)的排列和运动的模拟结果直接计算求和以实现宏观现象中的数值估算。可以直接模拟许多宏观现象,取得和实验相符合或可以比较的结

51、果,还可以提供微观结构、运动以及它们和体系宏观性质之间关系的极其明确的图象15。MD 以其不带近似、跟踪粒子轨迹、模拟结果准确16,而倍受研究者的关注,在物理、化学、材料、摩擦学等学科及纳米机械加工中得到广泛而成功的应用。3.4.2分子动力学的适用范围材料科学中计算模拟研究范围极为广泛,从埃量级的量子力学计算到连续介质层次的有限元或有限差分模型,可分为4个层次:电子、原子、显微组织和宏观层次(见图2-2)。MD主要是原子尺度上研究体系中与时间和温度有关的性质的模拟方法。图3-2 材料模型的层次划分最早将分子动力学方法用于材料研究中的是Vineyard于1960年探讨材料辐射损伤的动力学规律。模

52、拟结果给出了原子轨迹,这一工作使得过去对热力学性能的定性估计迈向对微观过程的定量研究。1964年Rahman用MD方法模拟液体氩,同时加进了周期性边界条件,结果他惊奇地发现可以用很少的粒子(864个)来反映真实系统的热力学性质。自此,凝聚态物质的分子动力学模拟成为可能,许多研究者纷纷投入这一研究工作。最初应用是基于偶函数,如Lennard2Jones势函数和Morse势函数,模型简单,运算量小,而得以在材料科学中广泛应用。但由于其未考虑到体积相关项,在研究材料的弹性系数性质和预言金属的结合能及空位形成能时,难以获得准确的结果18。EAM多体势15主要用于fcc型金属及其合金中,处理其结构、热力

53、学、表面、缺陷及液态金属等问题,也应用于hcp及bcc型金属及合金,以及半导体Si。一般,MD方法在中型机或微机上进行时,由于其内存和运算速度的限制,模拟研究只能限于5001000个原子的小系统。因而模拟结果虽然也能揭示一些微观结构的特征和规律,但与实际的大系统情况有较大差异。在并行处理系统上对更大量的原子系统进行模拟研究19,其结果必然会接近于实际,从而对生产实践将会更有实际指导意义。3.4.3分子动力学的应用（1）金属的液态结构在目前实验条件下,液态金属的结构及其变化尚很难精确测定。王鲁红等人20采用F2S型多体势描述了8种hcp型金属的液态微观结构并与实验相比较,模拟结果表明,Mg、Co

54、和Zn的势函数可以较好的描述其液态结构,Ti和Zr的势函数则不能;Be和Ru的势函数描述的液态结构较为合理,Hf则与一般液态结构特点不一致。李辉等人21采用EAM多体势模型,很好地描述了液态过渡金属Ni的结构变化特性。（2）高压相稳定性预测高压试验比较困难,但在近自由电子近似下,用赝势理论可以预测高压相稳定性。金朝晖22研究了高压下Mg的相变,在常压Mg以1×1012和6×1012K/s的速度冷却时获得hcp结构,冷却速度为1.2×1013K/s时,得到非晶结构。在高压(45GPa)下,Mg以8×1012和5×1013K/s的速度冷却后为bcc

55、结构,冷却速度为1×1014K/s时可得非晶结构。（3）热力学参数的计算刘洪波23利用EAM模型几等体积MD方法模拟计算Cu2Ni合金的Gibbs自由能,结果表明,计算的合金混合Gibbs自由能与实验值较一致。（4）薄膜形成过程薄膜研究是当今科学研究的热点之一。目前在很多薄膜制备方法中,都应用了低能离子轰击技术,在这些方法中,低能离子/表面相互作用在控制薄膜的微观结构方面起着重要作用。由于离子/表面相互作用发生在时间间隔小于10-12s内,因而特别适合于用MD方法对这一过程进行描述。（5）金属界面研究目前有关界面的MD研究的不多。罗旋24曾对Ag/N和Cu/Ni界面在弯曲状态下的力学

56、行为作了研究,结果给出的力学性能曲线与宏观规律相符合,这说明用MD方法研究界面问题是可行的。采用MD方法进行界面研究是十分有意义的工作,尤其是目前正急需解决的金属-陶瓷界面问题。（6）凝固及玻璃化转变的模拟由于实验室急冷技术所能达到的冷却速度一般<107K/s,对于许多合金和金属不能达到其玻璃转化所需的冷却速度,对玻璃态结构的研究受到了限制。用EAM框架下的有效配对势模拟了Pd的一级相变及玻璃转化过程25,得到的热力学参数可以同实验室进行比较,还得到了对应于结晶和玻璃转化的冷却速度。又如,常规的快冷技术冷却速度一般<107K/s,无法使Ni3Al非晶化。王鲁红26用计算方法知道在什么情况下能非晶化,弥补了实验的不足。另外,MD法还在固液相变27、固态相变28、表面能计算29、晶体生长30等几个方面得到应用。3.4.4分子动力学的最新进展最近,第一原理MD方法正在兴起,它对经典MD方法产生很大的挑战。它克服了经典MD方法中原子之间相互作用势不够精确的缺点。该方法是建立在密度函数理论和局域密度近似框架下的自洽伪势模拟方法。经典MD方法只能给出原子(离子)的位置以及晶体的结构,但第一原理M