周 佳， 耿 珺

(1.武漢大學(xué) 動(dòng)力與機(jī)械學(xué)院， 武漢 430072； 2.國(guó)核(北京)科學(xué)技術(shù)研究院有限公司， 北京 102209)

?

周 佳1， 耿 珺2*

(1.武漢大學(xué) 動(dòng)力與機(jī)械學(xué)院， 武漢 430072； 2.國(guó)核(北京)科學(xué)技術(shù)研究院有限公司， 北京 102209)

采用金屬銀的嵌入原子模型(EAM)，利用蒙特卡洛方法(MC)方法，在正則系綜(NVT)系綜下，計(jì)算了銀從1 700 K到2 300 K的飽和蒸氣壓，并和實(shí)驗(yàn)所測(cè)得的蒸汽壓計(jì)算公式進(jìn)行了對(duì)比.計(jì)算結(jié)果表明所有溫度下的模擬結(jié)果與實(shí)驗(yàn)測(cè)量的飽和蒸汽壓誤差均在30%以內(nèi)，驗(yàn)證了EAM勢(shì)能可以定性符合銀的飽和蒸氣壓，也證明了銀的EAM模型可以拓展到氣態(tài)的模擬.

飽和蒸汽壓; 蒙特卡洛模擬; 分子動(dòng)力學(xué); 相平衡

Daw和Baskes[2]針對(duì)有效介質(zhì)理論和準(zhǔn)原子理論存在的缺陷，提出了嵌入原子模型(Embedded Atom Method， EAM) .假設(shè)每一個(gè)原子皆為嵌入其他局部原子組成晶格中的客體，嵌入原子能即為嵌入該處的前后的能量之差，這個(gè)能量差是由平均電子密度決定的[3].最先使用這種方法來對(duì)純?cè)氐倪^渡金屬相關(guān)性質(zhì)進(jìn)行計(jì)算，包括結(jié)構(gòu)和熱力學(xué)性能，表面和晶界性能[2]，后EAM勢(shì)能推廣到了Ni-Al為代表的幾種fcc合金中[4-7]計(jì)算其熱力學(xué)性能以及拉伸性能，眾多結(jié)果表明EAM模型還能應(yīng)用到液態(tài)的計(jì)算中計(jì)算其原子輸運(yùn)性能[8].通過Johnson的截?cái)嘈拚齕9]，EAM模型也十分成功的描述了bcc與hcp金屬的熱力學(xué)和及拉伸性能[10-12]，另外還有基于該方法構(gòu)造的表面分析型模型計(jì)算[13].至今，基于EAM方法的氣液態(tài)研究較少，目前未見登載的使用該模擬方法計(jì)算飽和蒸汽壓的文獻(xiàn).本文通過編程完成了蒙特卡洛模擬(Monte Carlo Simulation，MC)，計(jì)算了飽和蒸汽壓.使用分子動(dòng)力學(xué)模擬(Molecular Dynamics，MD)驗(yàn)證了MC程序的正確性，并將計(jì)算出的飽和蒸汽壓結(jié)果與MB Panish[1]使用色譜法測(cè)量的Ag蒸汽壓進(jìn)行了研究比對(duì)，同時(shí)證明了EAM模型拓展到氣態(tài)的可行性，對(duì)EAM的研究領(lǐng)域進(jìn)行了拓展.

1 模型與方法

1.1 MC模擬

MC方法[14]是一種以概率統(tǒng)計(jì)理論為指導(dǎo)的數(shù)值計(jì)算方法.此次模擬使用C語言編寫了MC模擬程序，在編程過程中參考了Frenkel與Smit提出的模擬思路[15]，借助Metropolis抽樣算法[16]采樣.模擬計(jì)算輸出的固定溫度T下兩相壓強(qiáng)p與化學(xué)勢(shì)μ便是確定飽和蒸汽壓的依據(jù)，EAM勢(shì)能模型的具體公式為：

(1)

其中，Ei為原子i的能量，ρβ為原子i處原子j的β類型電子的密度，F(xiàn)為嵌入原子能，φαβ為原子i和原子j的兩體對(duì)勢(shì)，rij為i與j原子間的距離.勢(shì)能采用LAMMPS中提供的EAM勢(shì)能文件，參考Foiles[3]的EAM參數(shù).規(guī)定原子質(zhì)量數(shù)為107.87u，截?cái)嗑嚯xrc為5.5 ?，最大電子密度為0.25C/α0，將截?cái)嗑嚯x和電子密度均分為500段，列出了對(duì)應(yīng)的F，ρβ與φαβ.輸出的壓強(qiáng)采用Foiles提出的EAM壓強(qiáng)公式[17]，使用Widom插入粒子方法[18，19]，通過計(jì)算受驗(yàn)粒子與系統(tǒng)相互作用求化學(xué)勢(shì)：

(2)

式中，μex表示的是N個(gè)粒子的體系內(nèi)插入第N+1個(gè)粒子時(shí)增加的化學(xué)勢(shì)，n表示數(shù)密度，λ表示該溫度下的熱德布羅意波長(zhǎng)，ψ表示加入原子后增加的勢(shì)能.本次計(jì)算使用模擬正則系綜(NVT)，512個(gè)Ag原子，MC循環(huán)的總次數(shù)為4×106次，從初始狀態(tài)達(dá)到平衡態(tài)的次數(shù)為2×106次，采樣頻率為每10步一次.

1.2 MD模擬

飽和蒸汽壓的計(jì)算不能直接使用MD模擬，因?yàn)闅鈶B(tài)與液態(tài)的壓強(qiáng)在數(shù)量級(jí)上存在很大的差異，很難使用穩(wěn)定的氣液兩相數(shù)據(jù)通過擬合方法得到壓強(qiáng)-體積曲線.Bhattacharya[20]曾計(jì)算出了溫度在沸點(diǎn)以上時(shí)的Cu和Al的壓強(qiáng)-體積曲線，這種方法對(duì)其他溫度并不適用，本文僅使用MD方法來驗(yàn)證MC的計(jì)算結(jié)果確定編制程序的可靠性.使用軟件LAMMPS[21]，輸出壓強(qiáng)，總勢(shì)能與原子間的對(duì)勢(shì).以1 700K為例，如圖1所示為氣態(tài)的輸出，結(jié)果在密度較小時(shí)，與理想氣體十分接近.該狀態(tài)下的原子間隔較大作用力幾乎為0，近似為理想氣體模型.在比對(duì)化學(xué)勢(shì)判斷相平衡時(shí)，取理想氣體的化學(xué)勢(shì)與MC計(jì)算的化學(xué)勢(shì)比對(duì).圖2為液態(tài)的輸出，壓強(qiáng)與密度在穩(wěn)定液態(tài)密度附近呈成線性相關(guān).從兩幅圖的數(shù)據(jù)可以看出氣態(tài)與液態(tài)的壓強(qiáng)在數(shù)量級(jí)上有很大差異.

圖1 1 700 K下氣態(tài)MD的壓強(qiáng)與理想氣體Fig.1 Pressure of MD and ideal gas at 1 700 K

金屬Ag的晶格為面心立方(fcc)，模擬系統(tǒng)采用8×8×8的晶格盒子，共2 048個(gè)原子，使用周期邊界條件，通過設(shè)定不同晶格常數(shù)使各原子在設(shè)定密度下均勻分布.通過Nose-Hoover熱浴法[22-23]積分方法為Verlet形式[24]控制溫度.經(jīng)過測(cè)試，氣態(tài)時(shí)選定的步長(zhǎng)為20fs液態(tài)5fs.溫度從1 700K到2 300K，間隔100K，穩(wěn)定氣態(tài)密度與液態(tài)密度與隨著模擬溫度變化.每次計(jì)算弛豫8×104步后以總能量來判斷平衡，再對(duì)其進(jìn)行熱力學(xué)性質(zhì)的輸出，共運(yùn)行1.6×106步.

圖2 1 700 K下液態(tài)MD模擬的Ag壓強(qiáng)Fig.2 MD pressure of liquid Ag at 1 700 K

2 模擬與計(jì)算

使用MC程序計(jì)算了溫度1 700K至2 300K的粒子平均能量，壓強(qiáng)及化學(xué)勢(shì)，輸出了MC計(jì)算的接受比率.以2 300K輸出的這部分結(jié)果為例.2 300K下的MC程序計(jì)算得到液態(tài)部分的結(jié)果如表1所示，將化學(xué)勢(shì)的單位轉(zhuǎn)化為約化單位n×kB×T便于后面的蒸汽壓計(jì)算，每次的模擬的位移dr設(shè)置為0.5 ?.圖3～圖9為1700～2300 K液態(tài)下MC與MD壓強(qiáng)結(jié)果對(duì)比，可以看到：MC與MD的壓強(qiáng)比對(duì)結(jié)果均出現(xiàn)了200 bar左右的波動(dòng)，此量級(jí)的波動(dòng)是由于液態(tài)下原子間的距離近，壓強(qiáng)計(jì)算公式中維里項(xiàng)較大造成的.溫度較低時(shí)可以觀察到MC模擬所得到的結(jié)果較為穩(wěn)定波動(dòng)更小，線性關(guān)系明顯，高溫的計(jì)算結(jié)果的波動(dòng)則相對(duì)偏大.

表1 2300 K下MC模擬的輸出結(jié)果

圖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 模擬飽和蒸汽壓與實(shí)驗(yàn)飽和蒸汽壓Fig.10 Comparing SVP with experimental

另外MC程序的計(jì)算與MD的計(jì)算所取的原子個(gè)數(shù)不相同，MD模擬的原子個(gè)數(shù)為MC模擬的4倍，體系更大也更穩(wěn)定，所以MD結(jié)果的線性關(guān)系明顯.未來的研究應(yīng)將上述模型中的模擬體系與平均次數(shù)提高，減小誤差.通過MC輸出的數(shù)據(jù)可以得到液相的壓強(qiáng)與化學(xué)勢(shì)的線性方程，氣相部分如上文所示取理想氣體壓強(qiáng)與化學(xué)勢(shì).將氣液兩相的壓強(qiáng)與化學(xué)勢(shì)的進(jìn)行比對(duì)，化學(xué)勢(shì)與壓強(qiáng)均相等時(shí)便可得到此溫度下的飽和蒸汽壓，此次計(jì)算結(jié)果如表2所示.

表2 飽和蒸汽壓計(jì)算結(jié)果

如表2所示，所有溫度下的模擬結(jié)果與實(shí)驗(yàn)測(cè)量的飽和蒸汽壓誤差均在30%以內(nèi)，與實(shí)驗(yàn)值定性符合.溫度在1 700 K時(shí)的結(jié)果誤差最小為11.82%，最大的誤差則出現(xiàn)在1 800 K，為29.16%，1 900 ～2 100 K的溫度下的結(jié)果與實(shí)驗(yàn)值接近.1 900～2 300K模擬與實(shí)驗(yàn)值相比偏大，1 700 K與1 800 K則偏小.模擬計(jì)算隨著溫度的上升而增加，根據(jù)計(jì)算結(jié)果可得到飽和蒸汽壓與溫度的關(guān)系有：

其中， 溫度T的單位是K，p的單位是bar，其與溫度的關(guān)系如圖10所示，對(duì)比實(shí)驗(yàn)值，兩者較為接近.

3 結(jié)論

本次計(jì)算通過編程完成MC模擬，采用EAM模型得到了金屬Ag的飽和蒸汽壓，使用MD法對(duì)計(jì)算結(jié)果進(jìn)行驗(yàn)證，并且將計(jì)算結(jié)果與實(shí)驗(yàn)值[1]進(jìn)行了比對(duì).結(jié)果與實(shí)驗(yàn)測(cè)量誤差在30%以內(nèi)，證明了嵌入原子理論，合理成功地描述了金屬Ag的飽和蒸汽壓，該勢(shì)能模型可以從液態(tài)的模擬大量拓展到氣態(tài)模擬的研究與應(yīng)用.液態(tài)MC與MD模擬的飽和蒸汽壓對(duì)比的誤差，需要日后對(duì)模擬體系加以優(yōu)化.堿金屬的EAM勢(shì)能已經(jīng)十分完善[25]，今后將應(yīng)用這種方法研究對(duì)堿金屬和其他穩(wěn)定EAM勢(shì)能金屬的飽和蒸氣壓進(jìn)行驗(yàn)證.

[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] 張邦維. 嵌入原子方法理論及其在材料科學(xué)中的應(yīng)用:原子尺度材料設(shè)計(jì)理論[M]. 長(zhǎng)沙:湖南大學(xué)出版社， 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模型的發(fā)展現(xiàn)狀與展望[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ì)算[J]. 金屬學(xué)報(bào)， 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.

湖北省教育廳優(yōu)秀中青年人才科研項(xiàng)目(Q20114502)．

1000-1190(2015)03-0363-05

TG111.3

A

*通訊聯(lián)系人. E-mail: gengjun@snptc.com.cn.