题具体分析。对于一个复杂的系综，单个粒子所遵循的运动规律已经完全确定，多种类型的粒子间相互作用的力的形式也完全确定，这样对于一个粒子系综内部粒子的运动的描述具有明确的表达公式，但是其运动方程的解析却是不可能求得的。对于一小块样品，其内部粒子数也是很多的，
要想精确求出系综中粒子间相互作用和运动状态是不可能的，只有通过少量粒子的运动性质计算来求出宏观性质。由于粒子数过少，在构建的计算胞边缘的粒子较多，其受力状态与内部粒子不同，为了解决这一问题，分子动力学方法需要对计算胞设置适当的边界条件，主要有四种[11]：
(1)自由表面边界条件。计算胞中一般只包含几千个原子，对应10-20摩尔量级的材料。因此，自由表面边界条件适用于研究大的自由分子。
(2)刚性边界条件。计算胞的外围是一个具有晶体结构的原子位置被固定的“包层”。“包层”的厚度应大于原子间的作用距离，代表着与计算胞内原子有相互作用的宏观晶体。
(3)柔性边界条件。计算胞外围的原子在力作用下，可以有小的位移。研究扩展缺陷如位错运动时，常采用此种边界条件。
(4)周期性边界条件。一般用于模拟较大的系统，即在计算胞中，最左（或前、上）面的原子与最右（或后、下）面的原子之间有相互作用。
分子动力学方法的基本方程都是线性或非线性的二阶常微分方程，对此人们很早就提出了阿达姆斯(Adams)方法、欧拉(Euler)折线法、改进的欧拉法及其龙格-库塔(Runge-Kutta)方法等，并在实际中得到了广泛应用。就分子动力学而言，主要使用Verlet法和Gear法进行数值积分求解，具体求解方法可参考文献[12]。
所谓分子动力学模拟，是指对于原子核和电子所构成的多体系统，用计算机模拟原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在全部其它原子核核电子所提供的经验势场作用下按牛顿定律运动。
16
图3-1 分子动力学的基本原理流程
分子动力学模拟的基本原理：分子动力学将连续介质看成由N个原子或分子组成的粒子系统,各粒子之间的作用力可以通过量子力学势能函数求导得出,忽略量子效应后,运用经典牛顿力学建立系统粒子运动数学模型,通过数值求解得到粒子在相空间的运动轨迹,然后由统计物理学原理得出该系统相应的宏观动态、静态特性。图1所示是MD 模拟过程。MD具体的做法是计算机上求运动方程的数值解。通过适当的格式对方程进行近似,使之适于在计算机求数值解。从使用连续变数和微分算符的描述过渡到使用离散变数和有限差分算符的描述,显然会有误差,误差的阶数取决于具体的近似机制,即所用的算法。模拟首先是规定初始条件。为了确定起见,可令初始位置在格子的格点上,而初始速度则从波尔兹曼分布得出。一个按上述办法建立的系统不会具有所要的能量,而且,很可能这个状态并不对应于一个平衡态。为了推动系统到达平衡,需要一个趋衡阶段。可以通过增加或从系统中移走能量,对运动方程向前积分若干时间步,使系统弛豫到平衡态。接着是物理量的计算阶段,沿着系统在相空间中的轨道计算一切令人感兴趣的量。
模拟中,MD采用周期边界条件和最小镜像原理,可以大幅度减少计算工作量[14]。周期边界条件是将一定数量的粒子N集中在一定的容积V中,这个容积V称为原胞,原胞周围的部分可以看作是原胞的复制,它们称作镜像细胞。这些镜像细胞的尺寸和形状与原胞完全相同,并且每个镜像细胞所包含的N个粒子是原胞中粒子的镜像,原胞在各个方向周期复制便形成了宏观物质样本。这样只需根据原胞周围的边界条件计算原胞内粒子的运动,幅度减少了工作量。原子间作用势能模型的构造对于MD法的应用至关重要。最简单的偶势模型只考虑两体作用,而与其它原子无关,在模拟中运算量小。20
17
世纪80年代以来,各种经验或半经验的多体势模型迅速发展,特别是镶嵌原子(EAM)[17]既克服了偶势的缺陷,又不会使计算量太大。 3.3.2 原子间相互作用势
在分子动力学模拟中，能够正确的描述原子之间的相互作用势是至关重要的，它
决定了模拟结果的可靠性[13]，因此，分子动力学模拟的很大部分工作被放在了对原子间相互作用势的拟合上。将原子间的相互作用用一个简单的势函数表示，由势函数出发即可以得到物质的各种静态物理性质，如结合能、点阵常数和弹性模量等，这样就把问题大大简化了。通过原子间相互作用势对位移求导数就可以用来描述原子之间的相互作用力，根据前面提到的求解诸多粒子经典动力学方程完成分子动力学模拟，就可以获得物质的一些动态性质，如相稳定性等。基于量子力学的第一性原理计算原子间的相互作用势计算量大，对计算机硬件要求过高，因此在计算机模拟中，人们发展出了各种经验的和半经验的作用势进行使用。最先使用的是对势，后来在对势的基础上又发展出了几种多体势。势参数的拟合
在选取适当的势函数形式，并指定F(ρ),f(r)和φ(r)的函数形式后，我们需要针对自己所研究的金属系统拟合势参数，得到合适的多体势。首先给定势参数的初始值，其次利用该多体势导出该系统的物理性质，例如晶格常数，弹性模量等，并与实验值相比较，然后微调势参数，直到获得满意的结果。通常，还通过拟合的多体势计算能量曲线，与Rose方程进行比较以获得对势参数的进一步认可。对于一个由A、B 组成二元合金系统，一个完整的多体势应有三组势参数：即描述同种原子A-A、B-B 之间的作用势的势参数以及描述异种原子A-B之间的作用势（也称为交叉势）的势参数。因此，需要具有纯金属A、B以及AB所形成的化合物的一些物理性质数据用于势参数的确定。对于纯金属以及生成热为负的合金系统的化合物其数据一般可以通过工具书查找，但对于一些生成热为正的系统，由于化合物不能稳定存在，其物理性质数据难以得到。为了解决这个问题，研究组尝试利用第一性原理辅助构建多体势来解决这个问题，通过第一性原理计算得到合适的亚稳化合物及其物理性质数据，用于交叉势势参数的拟合，目前已经取得较好的效果。前面获得的多体势参数是通过对物质静态性能拟合得到的，为了确保多体势的准确性，还需通过分子动力学模拟验证一些动态性能，
18
如晶胞稳定性、相能量序、熔点等。如果利用构建的多体势模拟出来的动态性能与实际出入太大，则需要继续寻找新的势参数，直到静态性能和动态性能都与实际吻合较好，因此多体势的构建是一个极度耗时的实验校正（trial and error）过程。获得了一组优良的多体势参数后，就可以利用此多体势进行分子动力学模拟，得到自己所需要的结果。
3.4 分子动力学在材料科学中的应用
3.4.1分子动力学的概述
分子动力学(Molecular Dynamics,简MD) 用于计算以固体、液体、气体为模型的单个分子运动,它是探索各种现象本质和某些新规律的一种强有力的计算机模拟方法,具有沟通宏观特性与微观结构的作用,对于许多在理论分析和实验观察上难以理解的现象可以做出一定的解释[14]。MD方法不要求模型过分简化,可以基于分子(原子、离子)的排列和运动的模拟结果直接计算求和以实现宏观现象中的数值估算。可以直接模拟许多宏观现象,取得和实验相符合或可以比较的结果,还可以提供微观结构、运动以及它们和体系宏观性质之间关系的极其明确的图象[15]。MD 以其不带近似、跟踪粒子轨迹、模拟结果准确[16],而倍受研究者的关注,在物理、化学、材料、摩擦学等学科及纳米机械加工中得到广泛而成功的应用。 3.4.2分子动力学的适用范围
材料科学中计算模拟研究范围极为广泛,从埃量级的量子力学计算到连续介质层次的有限元或有限差分模型,可分为4个层次:电子、原子、显微组织和宏观层次(见图2-2)。MD主要是原子尺度上研究体系中与时间和温度有关的性质的模拟方法。
19
图3-2 材料模型的层次划分
最早将分子动力学方法用于材料研究中的是Vineyard于1960年探讨材料辐射损伤的动力学规律。模拟结果给出了原子轨迹,这一工作使得过去对热力学性能的定性估计迈向对微观过程的定量研究。1964年Rahman用MD方法模拟液体氩,同时加进了周期性边界条件,结果他惊奇地发现可以用很少的粒子(864个)来反映真实系统的热力学性质。自此,凝聚态物质的分子动力学模拟成为可能,许多研究者纷纷投入这一研究工作。最初应用是基于偶函数,如Lennard2Jones势函数和Morse势函数,模型简单,运算量小,而得以在材料科学中广泛应用。但由于其未考虑到体积相关项,在研究材料的弹性系数性质和预言金属的结合能及空位形成能时,难以获得准确的结果[18]。EAM多体势[15]主要用于fcc型金属及其合金中,处理其结构、热力学、表面、缺陷及液态金属等问题,也应用于hcp及bcc型金属及合金,以及半导体Si。一般,MD方法在中型机或微机上进行时,由于其内存和运算速度的限制,模拟研究只能限于500～1000个原子的小系统。因而模拟结果虽然也能揭示一些微观结构的特征和规律,但与实际的大系统情况有较大差异。在并行处理系统上对更大量的原子系统进行模拟研究[19],其结果必然会接近于实际,从而对生产实践将会更有实际指导意义。 3.4.3分子动力学的应用 （1）金属的液态结构
在目前实验条件下,液态金属的结构及其变化尚很难精确测定。王鲁红等人[20]采用F2S型多体势描述了8种hcp型金属的液态微观结构并与实验相比较,模拟结果表
20