明,Mg、Co和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结构,冷却速度为1×1014K/s时可得非晶结构。
（3）热力学参数的计算
刘洪波[23]利用EAM模型几等体积MD方法模拟计算Cu2Ni合金的Gibbs自由能,结果表明,计算的合金混合Gibbs自由能与实验值较一致。 （4）薄膜形成过程
薄膜研究是当今科学研究的热点之一。目前在很多薄膜制备方法中,都应用了低能离子轰击技术,在这些方法中,低能离子/表面相互作用在控制薄膜的微观结构方面起着重要作用。由于离子/表面相互作用发生在时间间隔小于10-12s内,因而特别适合于用MD方法对这一过程进行描述。 （5）金属界面研究
目前有关界面的MD研究的不多。罗旋[24]曾对Ag/N和Cu/Ni界面在弯曲状态下的力学行为作了研究,结果给出的力学性能曲线与宏观规律相符合,这说明用MD方法研究界面问题是可行的。采用MD方法进行界面研究是十分有意义的工作,尤其是目前正急需解决的金属-陶瓷界面问题。 （6）凝固及玻璃化转变的模拟
由于实验室急冷技术所能达到的冷却速度一般<107K/s,对于许多合金和金属不能达到其玻璃转化所需的冷却速度,对玻璃态结构的研究受到了限制。用EAM框架下的有效配对势模拟了Pd的一级相变及玻璃转化过程[25],得到的热力学参数可以同实验室进行比较,还得到了对应于结晶和玻璃转化的冷却速度。又如,常规的快冷技术冷却速度一般<107K/s,无法使Ni3Al非晶化。王鲁红[26]用计算方法知道在什么情况下能非晶
21
化,弥补了实验的不足。另外,MD法还在固液相变[27]、固态相变[28]、表面能计算[29]、晶体生长[30]等几个方面得到应用。 3.4.4分子动力学的最新进展
最近,第一原理MD方法正在兴起,它对经典MD方法产生很大的挑战。它克服了经典MD方法中原子之间相互作用势不够精确的缺点。该方法是建立在密度函数理论和局域密度近似框架下的自洽伪势模拟方法。经典MD方法只能给出原子(离子)的位置以及晶体的结构,但第一原理MD方法却能同时给出晶体结构及基态性质。第一原理分子动力学的应用较多的是分子族计算方法中的离散变分Xa法[15],其在材料科学中的应用主要是以下几个方面:合金化改善Ni3Al塑性的电子结构机制;材料硬度的第一原理标度;Co3Ti合金的环境脆性;Ni/Ni3Al界面结合的杂质效应;金属填充碳纳米管;合金元素在金属间化合物中的替代行为;抗热腐蚀镍基单晶合金的设计。
龚新高[31]用第一原理MD方法研究了高温下正常密度和高密度液体镓的结构和性质。正确地得到了高温下的液态镓的结构因子,随着压强(密度)的增大,观察到结构因子主峰右边肩膀的逐渐形成,这是以前采用其他方法不能得到的。
分子动力学模拟的历史演变随着算法的改进和计算机的发展而发展的。改进算法,提高模拟的速度;进行大规模分子动力学方法研究,扩大模拟系统中的粒子数,减小计算误差,以获得稳定的精确的模拟[16]。第一原理分子动力学模拟方法具有广阔的应用前景,目前还处于发展阶段,基础理论和实际应用上的困难尚未得到完全解决,有待量子理论学家、材料学家和工程学家的共同合作来解决。
3.5 Gulp软件在分子动力学模拟中的应用
在二十世纪后半期随着计算机的飞速发展，在计算机上借助一定的软件来模拟固体材料的计算材料学逐渐发展起来。最初人们只是尝试着通过模拟去计算和解释实验上得到的数据，后来逐步发展到通过探索材料未知的领域来指导实验以及自行设计新材料的更高层次。相应地，也就发展起来了许多用于模拟的软件。二十世纪90 年代初期，以英国人Gale 为首的工作小组推出了The General Utility LatticeProgram
22
(GULP) 软件[43、44]。GULP软件的主要基础是力场原理，它主要用来模拟固体材料的性质。GULP 软件不仅集合了其它一些分子动力学软件（如PARAPOCS、DMAREL、DLPOLY 等[45、46]）的功能，而且还开发出了一些新的用于模拟固体材料性能的独特程序。
GULP 软件有许多功能，如晶体性质的计算，拟合势函数，过渡态计算，缺陷计算，表面计算等。这里主要用到的是晶体性质的计算和拟合势函数的功能。 势函数的拟合
拟合原理：GULP 软件具有通过拟合能量曲面或材料性质来得到未知势参数的功能，总体上它提供了两种途径。一种方法是：通过拟合第一原理计算的结果来得到势参数；另一种方法是：通过拟合可信度高的实验值来得到势参数。判断拟合质量高低的一个重要参数是平方和值F。无论用的是那种方法，我们都要用平方和公式去检验。具体平方和F 表示为：
F?式中
allobservable???fcdc?fabs?2
fcdc和fabs分别是计算值和已知的用于拟合的数值；?是每一平方项的权重因子。
对每个拟合值来说，它的权重因子的选取依赖于许多因素，如所计算的数的重要性及可信赖度，比如说晶体结构就比弹性常数更重要更可信， 而且弹性常数是为对势函数特别敏感的量，一般独立的弹性常数的计算值与标准值符合得很好的情形很少。最后进一步检验拟合好坏的一个依据是被拟合势参数的个数一定要少于用于拟合的标准值的个数。
拟合能量表面：拟合能量表面是一种较烦琐的方法。先通过第一原理计算出每个晶体结构的能量然后根据能量表面拟合势参数。通常由势函数计算的能量值与第一原理计算值有一定的偏差，而GULP 中设置了能量修正值解决了这一问题。 普通拟合：普通拟合是经验性拟合中最普遍的一种放方法。具体方法是首先 确定一个能恰当地描述材料内部原子之间势模型，然后通过拟合已知的晶体材料的结构数据、力学性质、介电性质等的本征信息，得到合理的势参数。用于拟合的已知数据可以是实验数据也可以是第一原理的计算数据，但准确性一定要高。
同时拟合：普通拟合的优势是速度快，但拟合精度差。当拟合壳模型势的时候，我们
23
需要把作用在每个原子上的力最小化。一般情况下我们并不知道壳模型最外层壳的坐标，除非离子被固定在对称中心上。过去，人们把壳固定在核上面，这样做相当于人为地使系统的极化最小化了，导致壳模型的效果降低，因而拟合的效果不佳。对于这种情况GULP软件的同时拟合方法可以允许壳位置在拟合中变化。在拟合过程中有两种方法可以实现，一是在拟合过程中在每个点把壳的位置最小化，另一种是添加额外的壳位置坐标参数作为拟合值。一般拟合结构时，这两种方法是等价的。而当去拟合其它性质时虽然误差很小，但相对来说第二种方法的准确性低一些。
弛豫拟合：在传统的拟合过程中，因为每个人在具体拟合结构信息数据及很多实验数据时的判断标准不同，所以得到一个好的平方和值F 并不一定是个好的拟合。对此弛豫拟合的处理方法是在每个点都对结构进行优化并计算结构参数的变化。在这种方法中壳模型自然很精确，而且无需进行同时拟合。不过这种拟合有时计算量太大，所以可以先进行一般拟合再进行弛豫拟合。
3.6 运用GULP软件进行均方位移(MSD)数据分析
原子模拟中运用的程序是GULP软件和Origin8统计软件，在DOS界面下运行GULP软件对不同温度下的文件进行运算并优化，运算出的结果以RES格式的文件进行保存。之后运用统计软件对所有的文件进行统计，得出的原子均方位移（MSD）数据以TXT格式的文件进行保存，最后再运行Origin8软件，打开对应的数据文件，确定X轴和Y轴画出不同时间下原子的运动方式并进行分析。
24
第4章 模拟结果分析与讨论
虽然温度在0℃以下，人们也会感到冰是润滑的．Gábor A．Somorjai等实验发现[1,2]：即使在-30到-10℃，冰的表面也是处于高度无序的“类水面”（water-like surface，他们认为这就是人们感到冰润滑的原因）；随温度升高至融点，无序层从冰表面逐渐扩展到深处．由此，不妨猜想，金属表面在较高温度下（如接近熔点时），也可能有如此表现．
与纯金属和合金体内原子相比，金属表面原子有较多的自由度，表面原子所受束缚较体内原子要小，因此表面形成后其原子容易离开原始位置扩散到其它位置．这将影响表面的许多特性：表面吸附、表面腐蚀与氧化、表面催化与晶体生长和表面粗糙度等．在较高温度下，如电子器件的钎焊过程中，微电子工业某些薄膜的生长过程中
[3]
，或者化工行业金属参与的催化反应中[4]，对金属的表面结构和表面扩散机理的深
入理解就显得非常重要．
近来，人们用分子动力学模拟方法在原子尺度上对金属表面扩散和熔化做了大量研究[3,5-7]．对于表面熔化，这些工作侧重于研究金属表面作为一个整体的熔化过程和熔体的微观结构变迁，而对表面熔化前的原子结构和原子运动还知之甚少．为了验证上述猜想（金属表面低于熔点时可能存在类似液体的结构），也为了加深对金属表面熔化前原子结构和原子运动的理解，本文以HfO2为例，采用分子动力学方法和原子间嵌入势，在熔点以下特定的温度区间，对HfO2不同取向表面，特别是第一层原子的热运动现象进行模拟和分析，发现HfO2在特定温度下，第一层原子呈现的“类液”现象，不同取向表面的“类液”结构有不同的自扩散形式．
对于HfO2 (110)面，模拟构型为6×8×8结构，共384个原子，z方向上面6
25