蒙特卡洛模擬金屬Ag的飽和蒸汽壓

周 佳， 耿 珺

(1.武漢大學 動力與機械學院， 武漢 430072； 2.國核(北京)科學技術(shù)研究院有限公司， 北京 102209)

?

周 佳1， 耿 珺2*

(1.武漢大學 動力與機械學院， 武漢 430072； 2.國核(北京)科學技術(shù)研究院有限公司， 北京 102209)

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

飽和蒸汽壓; 蒙特卡洛模擬; 分子動力學; 相平衡

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

1 模型與方法

1.1 MC模擬

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

(1)

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

(2)

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

1.2 MD模擬

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

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

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

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

2 模擬與計算

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

表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 模擬飽和蒸汽壓與實驗飽和蒸汽壓Fig.10 Comparing SVP with experimental

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

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

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

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

3 結(jié)論

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

1000-1190(2015)03-0363-05

TG111.3

A

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