Au 快速凝固及晶体生长的分子动力学模拟

http://www.paper.edu.cn - 1 - Au 快速凝固及晶体生长的分子动力学模拟 李成祥， 孟庆元 1 ，杨立军 哈尔滨工业大学航天学院，哈尔滨 (150001) qymeng@public.hr.hl.cn 摘要： 采用 Embedded-Atom Method (EAM)作用势函数，利用分子动力学模拟了冷却速 度和压力对液态 Au 凝固过程的影响。结果表明在相同压力下，冷却速度快时，形成非晶体； 冷却速度慢时，形成晶体，冷却速度越慢，结晶温度越高，结晶进行得越充分，得到的结构 越稳定。在同一冷却速度下，压力对晶体的结晶温度影响较大，压力越大，结晶温度越高， 但结晶过程所用的时间不受压力影响。数值算例给出了不同冷却速度下的热力学及分子动力 学过程。最后用固液两层构型的方法模拟了 Au 的晶体生长过程。 关键词：分子动力学模拟，EAM 势函数，相变热力学，相变动力学，晶体生长 中图分类号：O342 文献标示码：A 1. 引言 人工晶体研究包括晶体结构性能、晶体生长、晶体材料应用等方面，而晶体生长研究是 人工晶体研究的基础。晶体生长是一种复杂的相变过程，迄今为止，人们尚难以为实际晶体 生长过程提出一个完备的、定量的理论，所以晶体生长理论研究力图从本质上揭示晶体生长 过程的基本规律，进而指导实际晶体制备技术。 通过对液态金属凝固过程的研究可以了解晶体的生长过程，在凝固过程中，无序体系的 结构很难由实验的方法来确定。近年来，随着计算机技术的迅速发展，分子动力学（Molecular Dynamics, MD）模拟成为研究新材料的一种重要方法，在晶体生长过程的模拟中也取得了 很大的成功[1，2]。因此，计算机模拟成为了解过冷液态和非晶态等结构信息的重要手段。描 述金属体系的镶嵌原子法（EAM）构造的多体势函数已被成功地应用于位错、界面等方面[3]。 本文采用 EAM 势函数，利用分子动力学方法模拟了金（Au）在快速凝固过程中的热力学参 数和动力学参数的变化，并且模拟了压力的变化对凝固过程的影响，最后用计算机图形处理 技术对晶体的生长过程也进行了描述。 2. 原子间相互作用势 原子间相互作用势能采用的是由 Finnis 和 Sinclair[4]根据镶嵌原子法提出的多体势，体 系的总势能为 ;)()( 2 1 ∑∑ += ≠ ij iji itot rρFE ϕ ∑= ≠ ij iji rf j )(ρ (1) 式中 iρ 是第 i个原子所在位置处的电子密度，即除去原子 i之外其它原子对此位置的电子 http://www.paper.edu.cn - 2 - 密度的贡献； )( iF ρ 是镶嵌能函数； )( ijj rf 是原子 i周围第 j个原子贡献的电子密度； ijr 是 i和 j两原子之间的距离； )( ijrϕ 是两体中心势，其中的镶嵌能函数，电子密度和两体势的 确定及参数的选择取自文献[5]。 3. 分子动力学数值模拟过程 3.1 Au 快速凝固过程的模拟 分子动力学模拟的对象是施加周期性边界条件的立方盒子中的 864个Au原子,盒子的边 长为 2.448nm。选择的模拟时间步长为 14101 −× s，外部压力为零。首先让体系在 1800K 下 弛豫 15000 个时间步长，得到平衡的液态，然后分别以 12105× ， 11105× ， 11101× ， 10108× 1−⋅ sK 的速度进行冷却。每隔 200K 记录一次体系构型，每一个构型都运行 1500 个 时间步长，从而得到偶分布函数数据；能量和体积数据每隔 5K 记录一次；平均平方位移数 据每隔 1K 记录一次。 3.2 压力对凝固过程的影响 结晶的过程中，在高压的条件下，经常能诱发成核，压力对成核和晶体生长过程有着重 要的影响。在上面的凝固过程中，选择降温速率为 110105 −⋅× sK ，分别模拟了 0， 2101× ， 4101× ， 5101× 个大气压的条件下的凝固过程。 3.3 晶体生长过程的模拟 模拟的对象是长方体盒子中的 768 个 Au 原子，x，y，z 方向分别为 2.448，6.528，1.632nm， 在三个方向施加周期性边界条件，y 方向上采用固-液-固交替扩展的模型[6]。原子的初始构 型为 y 方向上两头为晶态，各占原子总数的 1/4，是取自 400K 下弛豫 5000 个时间步长的构 型；中间的液态 Au 原子取自 1600K 下弛豫 5000 个时间步长的构型。模拟的时间步长取为 15101 −× s，将整个体系置于 600K 下运行 10000 个步长。将模拟过程中的原子坐标记录下来， 利用计算机图形处理技术，用动画的方式将每个原子显示出来，可以使人们直观地了解体系 结构的变化和 Au 的晶体生长过程。 4. 结构分析方法 4.1 原子的平均能量 体系的内能是由动能和势能两部分组成，而描述原子间相互作用的势函数是原子之间距 离的函数，所以体系的结构随温度发生变化的时候，其能量也将发生相应的变化，因此可以 通过计算单个原子的平均能量的变化来了解体系结构的变化。 4.2 偶分布函数 偶分布函数 g(r)定义为[7] 0)()( ρρ rrg = (2) 式中， )(rρ 表示距离某个原子 r处的粒子数密度； 0ρ 表示体系的平均密度。偶分布函数与 结构衍射实验得到的干涉函数互为 Fourier 变换，因此是理论和实验对照的基本依据，也是 描述液态和非晶态等无序体系的有效方法。 http://www.paper.edu.cn - 3 - 4.3 体系的体积 在模拟的过程中，要保持外界给定的压力不变，为了维持体系内压与外压的平衡，体系 的原胞体积要发生改变。体积随温度的降低而变化，当发生液固转变时，体积将发生急剧变 化，以此来判断体系结构的变化。 4.4 平均平方位移 平均平方位移 MSD（Mean Square Displacement）定义为[8] ))0()((1)( 1 22 ∑ = −= N i ii rtrN tr (3) 式中，N 是原子总数； )0(ir 是原子 i在 t=0 时刻的位移； )(tri 是原子 i在 t 时刻的位移。用 MSD 随时间的变化可以表征液态金属原子的扩散行为。 5. 模拟结果与分析 5.1 热力学结果 图 1 给出了在不同的降温速率条件下，单个原子的平均能量和体系的体积随温度的变化 曲线，p1，p2，p3，p4 分别表示 12105× ， 11105× ， 11101× ， 10108× 1−⋅ sK 的速度进行 冷却的过程。从图中可以看出，在冷却速度快的 p1 中，能量的变化是连续降低的，曲线没 有转折点，此过程中 Au 形成非晶体；而 p2，p3，p4 中，能量曲线分别在 450K，550K，620K 左右陡然下降，这表明结晶过程已发生。体系体积的变化规律和能量的变化规律相同。由 p2，p3，p4 可以看出，随着冷却速度减慢，发生结晶的温度在升高。在液态区域内，四个 过程中能量和体积重合，说明不受冷却速度的影响。还可以看出，随着降温速度的降低，曲 线在转折处下降的越陡，最终的能量越低，说明下降单位温度所放出的能量越多，结晶越充 分，体系越稳定。由于模拟过程中的降温速率非常快，体系状态的变化总是要滞后于温度的 变化，所以用分子动力学方法模拟凝固过程的时候，容易造成体系过冷，使得模拟过程中结 晶的形核温度低于纯 Au 结晶的实验温度。 图 1 不同冷却速度下单个原子平均能量(a)和体系体积(b)随温度的变化 Fig.1 The mean energy E(a) of single atom and the atom volume V(b) vs temperature T http://www.paper.edu.cn - 4 - 图 2 不同冷却速度下的偶分布函数 Fig.2 Pair correlation function at different cooling rate 图 3 相同冷却速度不同压力下单个原子平均能量和体系体积随温度的变化 Fig.3 The mean energy of single atom and the volume vs temperature under definite cooling rate 图 3 是在相同冷却速度不同压力下单个原子平均能量和体系体积随温度的变化，冷却速 度为 110105 −⋅× sK ， K1，K2，K3，K4 分别表示在施加外部压力为 0， 2101× ， 4101× ， 5101× 个大气压下的凝固过程。在 K1，K2，K3，K4 过程中，能量曲线分别在 540K，560K， 610K，850K 左右陡降，说明随着施加压力的增大，体系总能量降低，而结晶的温度逐渐升 高，从结晶温度的增加量可以看出，在纳米尺寸下，压力对晶体的结晶温度影响很大，压力 越大，越容易形核。而从体积变化曲线上可以看出，随着压力的增加，体系的初始体积由 17.5nm3 减小到 15.5nm3，而整个过程中的体积减少量由 2.5nm3 减小到 1.3nm3。能量曲线转 折处的斜率在结晶过程中没有变化，说明结晶过程所用的时间不受压力的影响。 5.2 动力学结果 平均平方位移（MSD）的变化能够反映凝固过程中原子的扩散行为，图 4 是 p1，p3， p4 三个过程中的 MSD 随温度的变化曲线。从图中可 以看出，p1 过程中的 MSD 变化较小，说明冷却速度 过快时，原子来不及发生大规模的移动就被“冻结” 在平衡位置附近振动，形成非晶体，p3，p4 过程的 MSD 变化就比较大，说明形成晶体的过程中原子作 了大规模的移动。由于使用了周期性边界条件，受 http://www.paper.edu.cn - 5 - 到原胞尺寸的限制，MSD 不能一直增大。 图 4 MSD 随温度的变化 Fig. 4 MSD vs the temperature 5.3 Au 晶体生长过程的分子动力学模拟 图 5 是 Au 晶体生长过程中 y-x 平面上的原子构型。(1)是初始构型，中间是液态 Au 原 子在 1600K 下弛豫 5000 个步长的构型，两端为晶态，各占原子总数的 1/4，是取自 400K 下 弛豫 5000 个步长的构型。将整个体系置于 600K 下运行 10000 个步长。在 t=0.5ps 时，固液 两种状态的原子相互作用，液态将固态部分原子融化，但随着时间的延长，液态部分又开始 慢慢结晶，由图(1)-(6)可以看到整个结晶过程的情况。 （1） t=0 ps (2) t=0.5 ps (3) t=2.5 ps (4) t=5.0 ps (5) t= 7.5 ps (6) t=9.0 ps 图 5 Au 晶体生长过程示意图(y-x 平面) Fig.5 representation of Au crystal growth (y-x plane) 6. 结论 1）采用 EAM 势函数，利用分子动力学方法模拟了 Au 的快速凝固过程，冷却速度快时 http://www.paper.edu.cn - 6 - 形成非晶体，冷却速度慢时形成晶体。冷却速度越慢，结晶温度越高，结晶进行得越充分， 得到的结构越稳定。 2）在纳米尺寸下，压力对晶体的结晶温度影响较大，压力越大，越容易形核，结晶温 度越高，但结晶过程所用的时间不受压力影响。 3）当冷却速度快时，平均平方位移变化较小，说明原子没有大规模扩散，形成的是非 晶体；当冷却速度慢时，平均平方位移变化较大，说明原子位置移动较大，形成晶体。 4）固液两层构型的方法能够很好地模拟 Au 的晶体生长过程。两端采用晶体有利于中 间液态部分的结晶形核，加快了晶体生长的过程。 参 考 文 献 [1] Broughton G Q and Gilmer G H. Molecular Dynamics Investigation of the Crystal-fluid III. Dynamical Properties of fcc Crystal-vapor Systems. J. Chem. Phys. ,1983,79:5119 [2] 李小平，韩其勇，刘洪波.金属学报，1995，8（31）：A356 [3] Foiles S M, Baskes M I, Daw M S. Phys Rev, 1986, B33: 7983 [4] G. J. Ackland，G. Tichy，V.Vitek and M. W. Finnis. Simple N-body Potentials for the Noble Metals and Nickel. Philosophical Magazine A. ,1987，Vol.56:735 [5] J. Mei, J.W. Davenport. Free-energy Calculations and the Melting Point of AI. Phys. Rev.., 1992，B46:21 [6] 陈柳，程兆年，唐鼎元. 晶体生长过程的分子动力学模拟研究. 人工晶体学报，1999，(28)：1096 [7] 李洪，谭志杰，张卫.外场导致的分形体的结构及偶分布函数.武汉大学学报，vol 46：353 [8] 王丽，边秀房，李辉.金属 Cu 液固转变及晶体生长的分子动力学模拟.物理化学学报，2000，(16)：825 The Molecular Dynamics Simulation of High-Speed Solidification and Crystal Growth of Gold (Au) Li Cheng-xiang Meng Qing-yuan Yang Li-jun Harbin Institute of Technology, Harbin (150001) Abstract The molecular dynamics simulation by using the embedded-atom method(EAM) was conducted to investigate how the solidification process of liquid gold Au is effected by the cooling rate and external pressure. The numerical simulation results show that the liquid Au is solidified into amorphous solid at high cooling rate and into crystal at relatively lower cooling rate under definite external pressure; The crystallization temperature becomes higher as the cooling rate goes down, in addition, the crystallization process is shown to be more sufficient and the molecular structure is more stable. Under a definite cooling rate, the pressure effects the crystallization temperature significantly. Higher the pressure, higher the crystallization http://www.paper.edu.cn - 7 - temperature. However, the crystallization time period is uneffected by the pressure. The numerical examples showed a complete thermodynamic and molecular dynamics process. Finally, the crystal growth process of liquid Au was simulated by using a two layer solid-liquid configuration model. Key words: The molecular dynamics simulation, EAM potential function, The phase transformation thermodynamics, The phase transformation dynamics, Crystal growth 作者简介：李成祥，男，1979 年生，博士研究生，主要研究方向为分子模拟。