时间步长对铜平衡熔点模拟计算结果的影响.pdf_大学文库

D01:10.13374j.isml00103x2006.06.010 第28卷第6期北京科技大学学报 Vol.28 Na 6 2006年6月 Journal of University of Science and Technology Beijing Jum.2006 时间步长对铜平衡熔点模拟计算结果的影响张远政倪晓东范韶蓉北京科技大学应用科学学院北京100083 摘要采用二分搜索等温等压分子动力学方法和EAM(embedded atom method)势，对单质铜的平衡熔点进行计算研究了不同时间步长对模拟计算结果的影响.结果表明：时间步长越小，熔点计算结果越接近实验值：当时间步长小于Q05s后，进一步减小时间步长，熔点计算结果在计算精度范围内没有变化关键词平衡熔点：时间步长：分子动力学：EAM势分类号TG146.1+1 材料熔化过程的计算模拟是分子动力学应用 (2)令的一个重要方面.这些研究的一个共同特点 T=I+Te (2) 就是采用升温来模拟材料的熔化过程，其熔点的 2 计算结果和升温过程有关2.理论上金属单质并让初始单胞在温度T下等温弛豫如果体的平衡熔点是材料的一个本征特性，其平衡熔化系不熔化则令T1=T:如果体系熔化则令T2 是在平衡熔点进行的等温过程，与材料加热升温 =T. 过程无关.为消除升温过程对平衡熔化计算的影 (3)反复进行以上过程直到T1和T2相差响，采用等温分子动力学驰豫和二分搜索法相结达到设定的计算精度要求，取T2为平衡熔点. 合的计算方法，研究铜的平衡熔化过程，同时考察计算单胞是一个施加周期性边界条件的6× 时间步长对平衡熔点计算模拟结果的影响. 6X6的完整立方晶体单胞，共864个Cu原子. 计算采用XMD分子动力学计算程序，Cu原子之 1 模拟计算方法间的相互作用利用Johnson的EAM势9描述，并首先构造一个含有N个原子的完整晶胞作采用Andersen方法7实现体系计算模拟中的等为初始计算单胞，在指定温度T下对初始计算单温等压控制，计算压强OPa.搜索过程初始温度胞进行等温等压弛豫直到体系达到平衡，然后通 T1为1000K,T2为2000K,每个温度下体系等过双体径向分布函数来判断体系是否熔化，温等压弛豫100000步，计算精度设为1K. 双体分布函数g(r)定义为： 2结果与讨论 g)品 (1) 图1~4给出了不同时间步长下，初始单胞在其中，p(r)表示距离任意一原子为r处的原子密不同温度等温弛豫后的双体分布函数曲线，图中度（统计平均值），P(0)是物质的平均原子密两相邻双体分布函数曲线沿竖轴相对平移一个单度9 位.由图1可以看出，当时间步长为5s时，温度二分搜索等温弛豫计算平衡熔点的过程为：低于或等于1688K的体系双体分布函数的峰值 (1)让初始计算单胞在一个较低的温度下弛明显，体现出了明显的晶体特征：温度等于或大于豫，如果体系不熔化则计该温度为T1:同样让初 1689K时双体分布函数第3峰和第4峰基本消始计算单胞在一个较高的温度下弛豫，如果体系失，体现出了明显的非晶特征，从而在1K的精度熔化则计该温度为T2· 范围内可以得到单质铜的平衡熔点为1689K.从图2~4则同样可以得到不同时间步长下铜的平收稿日期：2005-04-11修回日期：20050906 基金项目：国家自然科学基金资助项目(No.50472092) 衡熔点，结果列于表1. 作者简介：张远政(1980一，男，顾士研究生

时间步长对铜平衡熔点模拟计算结果的影响张远政倪晓东范韶蓉北京科技大学应用科学学院, 北京 100083 摘要采用二分搜索等温等压分子动力学方法和 EAM(embedded atom method)势, 对单质铜的平衡熔点进行计算, 研究了不同时间步长对模拟计算结果的影响.结果表明:时间步长越小, 熔点计算结果越接近实验值;当时间步长小于 0.05 fs 后, 进一步减小时间步长, 熔点计算结果在计算精度范围内没有变化. 关键词平衡熔点;时间步长;分子动力学;EAM 势分类号 TG 146.1 +1 收稿日期:2005 04 11 修回日期:2005 09 06 基金项目:国家自然科学基金资助项目(No .50472092) 作者简介:张远政(1980—), 男, 硕士研究生材料熔化过程的计算模拟是分子动力学应用的一个重要方面[ 1 4] .这些研究的一个共同特点就是采用升温来模拟材料的熔化过程 , 其熔点的计算结果和升温过程有关[ 2 4] .理论上金属单质的平衡熔点是材料的一个本征特性, 其平衡熔化是在平衡熔点进行的等温过程, 与材料加热升温过程无关 .为消除升温过程对平衡熔化计算的影响,采用等温分子动力学驰豫和二分搜索法相结合的计算方法, 研究铜的平衡熔化过程 ,同时考察时间步长对平衡熔点计算模拟结果的影响. 1 模拟计算方法首先构造一个含有 N 个原子的完整晶胞作为初始计算单胞 ,在指定温度 T 下对初始计算单胞进行等温等压弛豫直到体系达到平衡 ,然后通过双体径向分布函数来判断体系是否熔化. 双体分布函数 g(r)定义为: g(r)= ρ(r) ρ(0) (1) 其中 , ρ(r)表示距离任意一原子为 r 处的原子密度(统计平均值), ρ(0)是物质的平均原子密度 [ 5] . 二分搜索等温弛豫计算平衡熔点的过程为: (1)让初始计算单胞在一个较低的温度下弛豫,如果体系不熔化则计该温度为 T1 ;同样让初始计算单胞在一个较高的温度下弛豫 , 如果体系熔化则计该温度为 T2 . (2)令 T = T1 +T 2 2 (2) 并让初始单胞在温度 T 下等温弛豫, 如果体系不熔化, 则令 T1 =T ;如果体系熔化, 则令 T2 =T . (3)反复进行以上过程直到 T1 和 T2 相差达到设定的计算精度要求 ,取 T2 为平衡熔点. 计算单胞是一个施加周期性边界条件的 6 × 6×6 的完整立方晶体单胞, 共 864 个 Cu 原子. 计算采用 XMD 分子动力学计算程序,Cu 原子之间的相互作用利用 Johnson 的 EAM 势 [ 6] 描述,并采用 Andersen 方法[ 7] 实现体系计算模拟中的等温等压控制 , 计算压强 0 Pa.搜索过程初始温度 T 1 为 1 000 K , T2 为 2 000 K ,每个温度下体系等温等压弛豫 100 000 步 ,计算精度设为 1 K . 2 结果与讨论图 1 ～ 4 给出了不同时间步长下,初始单胞在不同温度等温弛豫后的双体分布函数曲线 ,图中两相邻双体分布函数曲线沿竖轴相对平移一个单位.由图 1 可以看出,当时间步长为 5 fs 时 ,温度低于或等于 1 688 K 的体系双体分布函数的峰值明显 ,体现出了明显的晶体特征;温度等于或大于 1 689 K 时双体分布函数第 3 峰和第 4 峰基本消失,体现出了明显的非晶特征 ,从而在 1 K 的精度范围内可以得到单质铜的平衡熔点为 1689K .从图 2 ～ 4 则同样可以得到不同时间步长下铜的平衡熔点,结果列于表 1 . 第 28 卷第 6 期 2006 年 6 月北京科技大学学报 Journal of University of Science and Technology Beijing Vol .28 No.6 Jun.2006 DOI :10.13374/j .issn1001 -053x.2006.06.010

。552· 北京科技大学学报 2006年第6期 I405K 689 1397K 688K 1396K 680K 1390K 0.2 0.40.60.81.012 0.2 0.40.60.810 1.2 rnm r/nm 图1时间步长为5时体系在不同温度下的双体分布函数图4时间步长为Q.005fs时体系在不同温度下的双体分布 Fig.I Pair distribution function g(r)of Cu at different tem- 函数 peratures (Time step is 5fs) Fig.4 Pair distribution function g(r)of Cu at different tem perature (Time step is0005 fs) 由表1可以看出各时间步长下的模拟结果均高于实验值1356K9,这可能是由于采用近似的 565 相互作用势，模拟体系的粒子个数有限，以及模型 553 未考虑实际晶体中存在的表面、空位、位错等缺陷 1554K 引起的. 1545K 二分搜索等温弛豫方法得到的铜的熔点范围小于等于1K,而升温过程计算，铜的熔点范围大 02 0.4 0.60.8 1.0 1.2 于10K,文献2,4中熔化温度的计算值分别为 r/nm 1520和1600K,相对误差为10%及18%，本工作图2时间步长为Q5s时体系在不同温度下的双体分布函计算值为1396K,相对误差仅为3%.由表1还数可以看出时间步长对模拟计算的结果有直接的影 Fig.2 Pair distribution function g(r)of Cu at different tem- peratures (Time step is 0.5fs) 响，时间步长越小，模拟结果越接近实验值，随着时间步长的进一步减小模拟结果趋于恒定，这主要是以下两方面的原因：(1)时间步长越小，计算精度越高，模拟的结果就应该越准确，当时间步长小到一定程度后，时间步长变化对计算结果精度 405 的影响趋于稳定；(2)从模拟对应的实际物理过程 396 考虑，时间步长越小，重新标度温度的时间间隔越 395k 短，这表明体系与外界热交换程度越好，更接近平 385K 衡熔化过程，时间步长小到一定程度后，体系与外 03 0.40.6081.0 界的热交换足够良好，时间步长变化对热交换速 rinm 度的影响也趋于稳定，图3时间步长为05s时体系在不同温度下的双体分布函 3结论数 Fig.3 Pair distribution function g(r)of Cu at different tem (1)时间步长越小，二分搜索等温过程计算 perature (Time step is 0.05 fs) 的熔点值越接近实验值. 表1铜的平衡熔点计算结果 (2)当时间步长小于0.05fs后，进一步减小 Table 1 Cal culated results of the equilibrium melting point of Cu 时间步长，熔点计算结果在计算精度范围内没有时间步长/s 熔点温度/K 时间步长s熔点温度/K 变化. 5X10T5 1689 5X10-17 1396 参考文献 5X1016 1555 5X10-18 1397 [I]Shu D j Sun D Y.A molecular-dynamics study of the

图 1 时间步长为 5 fs时体系在不同温度下的双体分布函数 Fig.1 Pair distribution function g(r)of Cu at different temperatures(Time step is 5 fs) 图 2 时间步长为 0.5 fs 时体系在不同温度下的双体分布函数 Fig.2 Pair distribution function g(r)of Cu at different temperatures(Time step is 0.5 fs) 图 3 时间步长为 0.05 fs时体系在不同温度下的双体分布函数 Fig.3 Pair distribution function g(r)of Cu at different temperature (Time step is 0.05 fs) 表 1 铜的平衡熔点计算结果 Table 1 Cal culated results of the equilibrium melting point of Cu 时间步长/s 熔点温度/ K 时间步长/ s 熔点温度/ K 5×10 -15 1 689 5×10 -17 1 396 5×10 -16 1 555 5×10 -18 1 397 图 4 时间步长为 0.005 fs 时体系在不同温度下的双体分布函数 Fig.4 Pair distribution function g(r)of Cu at different temperature (Time step is 0.005 fs) 由表 1 可以看出各时间步长下的模拟结果均高于实验值 1 356 K [ 8] ,这可能是由于采用近似的相互作用势,模拟体系的粒子个数有限 ,以及模型未考虑实际晶体中存在的表面、空位、位错等缺陷引起的[ 9] . 二分搜索等温弛豫方法得到的铜的熔点范围小于等于 1 K ,而升温过程计算 , 铜的熔点范围大于10 K [ 4] ,文献[ 2 , 4] 中熔化温度的计算值分别为 1 520和 1 600 K ,相对误差为 10 %及 18 %,本工作计算值为 1 396 K , 相对误差仅为 3 %.由表 1 还可以看出时间步长对模拟计算的结果有直接的影响,时间步长越小 ,模拟结果越接近实验值 ,随着时间步长的进一步减小模拟结果趋于恒定 .这主要是以下两方面的原因:(1)时间步长越小 ,计算精度越高 ,模拟的结果就应该越准确,当时间步长小到一定程度后, 时间步长变化对计算结果精度的影响趋于稳定;(2)从模拟对应的实际物理过程考虑 ,时间步长越小 ,重新标度温度的时间间隔越短,这表明体系与外界热交换程度越好 ,更接近平衡熔化过程,时间步长小到一定程度后 ,体系与外界的热交换足够良好 , 时间步长变化对热交换速度的影响也趋于稳定 . 3 结论 (1)时间步长越小 , 二分搜索等温过程计算的熔点值越接近实验值. (2)当时间步长小于 0.05 fs 后 ,进一步减小时间步长, 熔点计算结果在计算精度范围内没有变化 . 参考文献 [ 1] S hu D j, Sun D Y.A molecular-dynamics study of the · 552 · 北京科技大学学报 2006 年第 6 期