时间步长对铜平衡熔点模拟计算结果的影响

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 tem￾peratures(Time step is 5 fs) 图 2 时间步长为 0.5 fs 时体系在不同温度下的双体分布函 数 Fig.2 Pair distribution function g(r)of Cu at different tem￾peratures(Time step is 0.5 fs) 图 3 时间步长为 0.05 fs时体系在不同温度下的双体分布函 数 Fig.3 Pair distribution function g(r)of Cu at different tem￾perature (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 tem￾perature (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 期

Vol.28 No.6 张远政等：时间步长对铜平衡熔点模拟计算结果的影响 ·553。 arisotropic surfacemelting properties of AL (110).Surf Sci to applications.北京：化学工业出版社.2002 1999,441:206 [6 Johnson R A.Analytic nearest-neighbor model for foc metals. [2]Liu X,Meng C G.Melting and Superheating of Ag at High Phys Rev B 1988.37:3924 Heating Rate物理化学学报.20042D.280 [7]Andersen H C.Molecuhr dymamics simulations at constant [3]Zhang T.Zhang X R.Molecular dymamics simulation of the pressure and/or temperature.J Chem Phys 1980.72:2384 heating and melting processes of meal Cu.物理化学学报。 [8 Kreith F.The CRC Handbook of Themal Engineering. 2003,19:709 Bedin:CRC Press.2000 [4]Wang L Bian X F.Molecular dynamics simulat ion of metal [9 Lutsko J F.Wolf D Phillpot S R.Moleculr-dy namics study Cu during melting and crystallizing prcess..化学物理学报. of lattice-defect-nudeaed melting in metals using an embed 200013:544 ded-atom-method potential.Phys Rev B 1989.B40:2841 [5]Smit F.Understanding molecular simulation-from algorithms Influence of time step on the calculated melting point of Cu ZHANG Yuanzheng,NI Xiaodong,FAN Shaorong Appled Science School University of Science and Technology Beijing.Beijing 100083.China ABSTRACT Using temperature optimization process and embedded atom method (EAM),the equilibri- um melting point of Cu was calculated under the constant-temperature and pressure molecular dynamics. The influence of time step on the calculated results was also studied.The simulat ion shows that the shorter the time step is,the smaller the difference between the calculated melting point and the experimental is. When the time step is below 0 05 fs,the decreasing of time step has no effect on the calculated melting point in the range of calculation precision. KEY WORDS melting point;time step;molecular dynamics;EAM potential

anisotropic surface-melting properties of AL(110).Surf Sci, 1999 , 441:206 [ 2] Liu X , Meng C G .Melting and Superheating of Ag at High Heating Rat e.物理化学学报, 2004 , 20:280 [ 3] Zhang T , Zhang X R.Molecular dynamics simulation of the heating and melting processes of met al Cu.物理化学学报, 2003 , 19:709 [ 4] Wang L, Bian X F .Molecular dynamics simulation of metal Cu during melting and cryst allizing process.化学物理学报, 2000 , 13:544 [ 5] Smit F .Underst anding molecular simulation-from algorithms to applications.北京:化学工业出版社, 2002 [ 6] Johnson R A .Analytic nearest-neighbor model f or f cc metals. Phys Rev B, 1988 , 37:3924 [ 7] Andersen H C .Molecular dynamics simulations at const ant pressure and/ or t emperature.J Chem Phys, 1980 , 72:2384 [ 8] Kreith F .The CRC Handbook of Thermal Engineering . Berlin:CRC Press, 2000 [ 9] Lutsko J F, Wolf D, Phillpot S R.Molecular-dynamics study of lattice-def ect-nucleat ed melting in metals using an embed￾ded-at om-method potential.Phys Rev B, 1989 , B40:2841 Influence of time step on the calculated melting point of Cu ZHANG Y uanzheng , NI X iaodong , FAN Shaorong Applied S cience School, University of Science and T echnology Beijing , Beijing 100083 , C hina ABSTRACT Using temperature optimization process and embedded atom method (EAM ), the equilibri￾um melting point of Cu w as calculated under the constant-temperature and pressure molecular dy namics. The influence of time step on the calculated results was also studied .The simulation show s that the shorter the time step is, the smaller the difference between the calculated melting point and the experimental is. When the time step is below 0.05 fs, the decreasing of time step has no effect on the calculated melting point in the range of calculation precision . KEY WORDS melting point ;time step ;molecular dynamics ;EAM potential Vol.28 No.6 张远政等:时间步长对铜平衡熔点模拟计算结果的影响 · 553 ·