蒙特卡洛模拟金属Ag的饱和蒸汽压

周 佳， 耿 珺

(1.武汉大学 动力与机械学院， 武汉 430072； 2.国核(北京)科学技术研究院有限公司， 北京 102209)



周 佳1， 耿 珺2*

(1.武汉大学 动力与机械学院， 武汉 430072； 2.国核(北京)科学技术研究院有限公司， 北京 102209)

采用金属银的嵌入原子模型(EAM)，利用蒙特卡洛方法(MC)方法，在正则系综(NVT)系综下，计算了银从1 700 K到2 300 K的饱和蒸气压，并和实验所测得的蒸汽压计算公式进行了对比.计算结果表明所有温度下的模拟结果与实验测量的饱和蒸汽压误差均在30%以内，验证了EAM势能可以定性符合银的饱和蒸气压，也证明了银的EAM模型可以拓展到气态的模拟.

饱和蒸汽压; 蒙特卡洛模拟; 分子动力学; 相平衡

Daw和Baskes[2]针对有效介质理论和准原子理论存在的缺陷，提出了嵌入原子模型(Embedded Atom Method， EAM) .假设每一个原子皆为嵌入其他局部原子组成晶格中的客体，嵌入原子能即为嵌入该处的前后的能量之差，这个能量差是由平均电子密度决定的[3].最先使用这种方法来对纯元素的过渡金属相关性质进行计算，包括结构和热力学性能，表面和晶界性能[2]，后EAM势能推广到了Ni-Al为代表的几种fcc合金中[4-7]计算其热力学性能以及拉伸性能，众多结果表明EAM模型还能应用到液态的计算中计算其原子输运性能[8].通过Johnson的截断修正[9]，EAM模型也十分成功的描述了bcc与hcp金属的热力学和及拉伸性能[10-12]，另外还有基于该方法构造的表面分析型模型计算[13].至今，基于EAM方法的气液态研究较少，目前未见登载的使用该模拟方法计算饱和蒸汽压的文献.本文通过编程完成了蒙特卡洛模拟(Monte Carlo Simulation，MC)，计算了饱和蒸汽压.使用分子动力学模拟(Molecular Dynamics，MD)验证了MC程序的正确性，并将计算出的饱和蒸汽压结果与MB Panish[1]使用色谱法测量的Ag蒸汽压进行了研究比对，同时证明了EAM模型拓展到气态的可行性，对EAM的研究领域进行了拓展.

1 模型与方法

1.1 MC模拟

MC方法[14]是一种以概率统计理论为指导的数值计算方法.此次模拟使用C语言编写了MC模拟程序，在编程过程中参考了Frenkel与Smit提出的模拟思路[15]，借助Metropolis抽样算法[16]采样.模拟计算输出的固定温度T下两相压强p与化学势μ便是确定饱和蒸汽压的依据，EAM势能模型的具体公式为：

(1)

其中，Ei为原子i的能量，ρβ为原子i处原子j的β类型电子的密度，F为嵌入原子能，φαβ为原子i和原子j的两体对势，rij为i与j原子间的距离.势能采用LAMMPS中提供的EAM势能文件，参考Foiles[3]的EAM参数.规定原子质量数为107.87u，截断距离rc为5.5 Å，最大电子密度为0.25C/α0，将截断距离和电子密度均分为500段，列出了对应的F，ρβ与φαβ.输出的压强采用Foiles提出的EAM压强公式[17]，使用Widom插入粒子方法[18，19]，通过计算受验粒子与系统相互作用求化学势：

(2)

式中，μex表示的是N个粒子的体系内插入第N+1个粒子时增加的化学势，n表示数密度，λ表示该温度下的热德布罗意波长，ψ表示加入原子后增加的势能.本次计算使用模拟正则系综(NVT)，512个Ag原子，MC循环的总次数为4×106次，从初始状态达到平衡态的次数为2×106次，采样频率为每10步一次.

1.2 MD模拟

饱和蒸汽压的计算不能直接使用MD模拟，因为气态与液态的压强在数量级上存在很大的差异，很难使用稳定的气液两相数据通过拟合方法得到压强-体积曲线.Bhattacharya[20]曾计算出了温度在沸点以上时的Cu和Al的压强-体积曲线，这种方法对其他温度并不适用，本文仅使用MD方法来验证MC的计算结果确定编制程序的可靠性.使用软件LAMMPS[21]，输出压强，总势能与原子间的对势.以1 700K为例，如图1所示为气态的输出，结果在密度较小时，与理想气体十分接近.该状态下的原子间隔较大作用力几乎为0，近似为理想气体模型.在比对化学势判断相平衡时，取理想气体的化学势与MC计算的化学势比对.图2为液态的输出，压强与密度在稳定液态密度附近呈成线性相关.从两幅图的数据可以看出气态与液态的压强在数量级上有很大差异.

图1 1 700 K下气态MD的压强与理想气体Fig.1 Pressure of MD and ideal gas at 1 700 K

金属Ag的晶格为面心立方(fcc)，模拟系统采用8×8×8的晶格盒子，共2 048个原子，使用周期边界条件，通过设定不同晶格常数使各原子在设定密度下均匀分布.通过Nose-Hoover热浴法[22-23]积分方法为Verlet形式[24]控制温度.经过测试，气态时选定的步长为20fs液态5fs.温度从1 700K到2 300K，间隔100K，稳定气态密度与液态密度与随着模拟温度变化.每次计算弛豫8×104步后以总能量来判断平衡，再对其进行热力学性质的输出，共运行1.6×106步.

图2 1 700 K下液态MD模拟的Ag压强Fig.2 MD pressure of liquid Ag at 1 700 K

2 模拟与计算

使用MC程序计算了温度1 700K至2 300K的粒子平均能量，压强及化学势，输出了MC计算的接受比率.以2 300K输出的这部分结果为例.2 300K下的MC程序计算得到液态部分的结果如表1所示，将化学势的单位转化为约化单位n×kB×T便于后面的蒸汽压计算，每次的模拟的位移dr设置为0.5 Å.图3～图9为1700～2300 K液态下MC与MD压强结果对比，可以看到：MC与MD的压强比对结果均出现了200 bar左右的波动，此量级的波动是由于液态下原子间的距离近，压强计算公式中维里项较大造成的.温度较低时可以观察到MC模拟所得到的结果较为稳定波动更小，线性关系明显，高温的计算结果的波动则相对偏大.

表1 2300 K下MC模拟的输出结果

图3 1 700 K下的MC与MD输出Fig.3 MC and MC Ag results at 1 700 K

图4 1 800 K下的MC与MD输出Fig.4 MC and MC Ag results at 1 800 K

图5 1 900 K下的MC与MD输出Fig.5 MC and MC Ag results at 1 900 K

图6 2 000 K下MC与MD输出Fig.6 MC and MC Ag results at 2 000 K

图7 2 100 K下MC与MD输出Fig.7 MC and MC Ag results at 2 100 K

图8 2 200 K下的MC与MD输出Fig.8 MC and MC Ag results at 2 200 K

图9 2 300 K下的MC与MD输出Fig.9 MC and MC Ag results at 2 300 K

图10 模拟饱和蒸汽压与实验饱和蒸汽压Fig.10 Comparing SVP with experimental

另外MC程序的计算与MD的计算所取的原子个数不相同，MD模拟的原子个数为MC模拟的4倍，体系更大也更稳定，所以MD结果的线性关系明显.未来的研究应将上述模型中的模拟体系与平均次数提高，减小误差.通过MC输出的数据可以得到液相的压强与化学势的线性方程，气相部分如上文所示取理想气体压强与化学势.将气液两相的压强与化学势的进行比对，化学势与压强均相等时便可得到此温度下的饱和蒸汽压，此次计算结果如表2所示.

表2 饱和蒸汽压计算结果

如表2所示，所有温度下的模拟结果与实验测量的饱和蒸汽压误差均在30%以内，与实验值定性符合.温度在1 700 K时的结果误差最小为11.82%，最大的误差则出现在1 800 K，为29.16%，1 900 ～2 100 K的温度下的结果与实验值接近.1 900～2 300K模拟与实验值相比偏大，1 700 K与1 800 K则偏小.模拟计算随着温度的上升而增加，根据计算结果可得到饱和蒸汽压与温度的关系有：

其中， 温度T的单位是K，p的单位是bar，其与温度的关系如图10所示，对比实验值，两者较为接近.

3 结论

本次计算通过编程完成MC模拟，采用EAM模型得到了金属Ag的饱和蒸汽压，使用MD法对计算结果进行验证，并且将计算结果与实验值[1]进行了比对.结果与实验测量误差在30%以内，证明了嵌入原子理论，合理成功地描述了金属Ag的饱和蒸汽压，该势能模型可以从液态的模拟大量拓展到气态模拟的研究与应用.液态MC与MD模拟的饱和蒸汽压对比的误差，需要日后对模拟体系加以优化.碱金属的EAM势能已经十分完善[25]，今后将应用这种方法研究对碱金属和其他稳定EAM势能金属的饱和蒸气压进行验证.

[1] Panish B M. Vapor pressure of silver[J]. J Chem Eng Data， 1961(4):592-594．

[2] Baskes D . Embedded-atom method: Derivation and application to impurities， surfaces， and other defects in metals[J]. Physical Review B， 1984， 29(12):6443-6453．

[3] 张邦维. 嵌入原子方法理论及其在材料科学中的应用:原子尺度材料设计理论[M]. 长沙:湖南大学出版社， 2003.

[4] Foiles S M， Baskes M I， Daw M S. Embedded-atom-method functions for the fcc metals Cu， Ag， Au， Ni， Pd， Pt， and their alloys.[J]. Phys Rev B， 1986， 33(12):7983-7991．

[5] Foiles S M. Calculation of the surface segregation of Ni-Cu alloys with the use of the embedded-atom method[J]. Phys Rev B， 1985， 32(12):7685-7693．

[6] Johnson A R. Alloy models with the embedded-atom method.[J]. Physical Review B， 1989(17):12554-12559．

[7] Mei J， Davenport J W， Fernando G W. Analytic embedded-atom potentials for fcc metals: Application to liquid and solid copper[J]. Physical Review B， 1991， 43(6): 4653-4658．

[8] Asta M， Morgan D J， Hoyt J， et al. Embedded-atom-method study of structural， thermodynamic， and atomic-transport properties of liquid Ni-Al alloys[J]. Physical Review. B， 1999(22): 14271-14281．

[9] Johnson R A. Phase stability of fcc alloys with the embedded-atom method. [J]. Phys Rev B， 1990(41): 9717-9720．

[10] 胡望宇， 张邦维， 黄伯云. 分析型EAM模型的发展现状与展望[J]. 稀有金属材料与工程， 1999(1):1-4．

[11] Yifang O， Bangwei Z， Shuzhi L， et al. A simple analytical EAM model for bcc metals including Cr and its application[J]. Zeitschrift für Physik B， 1996， 101(2):161-168．

[12] Wu Y， Hu W， Sun L. Elastic constants and thermodynamic properties of Mg-Pr， Mg-Dy， Mg-Y intermetallics with atomistic simulations[J]. Journal of Physics D: Applied Physics， 2007， 40(23):7584-7592．

[13] 邓辉球， 胡望宇， 舒小林， 等. Pt-Rh二元合金系表面偏聚的分析型EAM模型计算[J]. 金属学报， 2001， 37(5): 467-471．

[14] Wasserstein R L， Wasserstein R L. Monte-Carlo: concepts， algorithms， and applications[J]. Technometrics， 1997， 39(3):338-350．

[15] Frenkel D， Smit B. Understanding molecular simulation[J]. Computers in Physics， 1997， 11(4):23-58

[16] Metropolis， Rosenbluth ， Rosenbluth ， et al. Equation of state calculations by fast computing machines[J]. The Journal of Chemical Physics， 1953， 21(6):1087-1092．

[17] Foiles M S. Application of the embedded-atom method to liquid transition metals.[J]. Phys Rev B， 1985， 32(6): 3409-3415．

[18] Widom B. Some topics in the theory of fluids[J]. The Journal of Chemical Physics， 1963， 39(11): 2808-2812．

[19] Binder K. Applications of Monte-Carlo methods to statistical physics[J]. Reports on Progress in Physics， 1997， 60(5):487-559．

[20] Bhattacharya C. Liquid-vapor phase diagram of metals using EAM potential[J]. AIP Conference Proceedings， 2013(1):52-53．

[21] Plimpton S. Fast parallel algorithms for short-range molecular dynamics[J]. Journal of Computational Physics, 1995, 117(1): 1-19.

[22] Nosé S. A unified formulation of the constant temperature molecular dynamics methods[J]. The Journal of Chemical Physics， 1984， 81(1):511-519．

[23] Hoover W G. Canonical dynamics: Equilibrium phase-space distributions[J]. Physical Review A， 1985， 31(3):1695-1697．

[24] Levesque D， Verlet L. Computer “experiments” on classical fluids (III). Time-dependent self-correlation functions[J]. Physical Review A， 1970 (6): 2514-2528．

[25] Belashchenko K D. Embedded atom method potentials for alkali metals[J]. Inorganic Materials， 2012， 48(1):79-86．

Monte-Carlo simulation of saturated vapor pressure of silver

ZHOU Jia1， GENG Jun2

(1.School of Power and Mechanical Engineering， Wuhan University， Wuhan 430072;2.State Nuclear Power Research Institute， Beijing 102209)

By means of embedded atom method(EAM) combined with the Monte Carlo and molecular dynamics simulation in the canonical ensemble (NVT)， saturated vapor pressure of silver was calculated from 1 700 K to 2 300 K， and the simulation results and vapor pressure calculation formula were compared. Results showed that the minimum error presented of 1 700 K. The error between simulation and the experimental results of all the temperature were within 30%. The EAM potential was verified to be able to accurately describe the liquid silver. saturated vapor pressure， which also deminstrated that the EAM model was able to be extended to the gaseous simulation from liquid simulation．

saturated vapor pressure; Monte-Carlo simulation; molecular dynamics; phase equilibrium

2015-01-20.

湖北省教育厅优秀中青年人才科研项目(Q20114502)．

1000-1190(2015)03-0363-05

TG111.3

A

*通讯联系人. E-mail: gengjun@snptc.com.cn.