搜索

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

面心立方Ce同构相变的分子动力学模拟

第伍旻杰 胡晓棉

引用本文:
Citation:

面心立方Ce同构相变的分子动力学模拟

第伍旻杰, 胡晓棉

Isostructural phase transition of fcc Ce: Molecular dynamics simulations

Diwu Min-Jie, Hu Xiao-Mian
PDF
HTML
导出引用
  • 基于嵌入原子法, 本文给出了一个金属Ce原子间的相互作用势. 利用该势分别计算了γ-Ce和α-Ce的晶格常数、结合能、弹性常数, 计算结果与实验或第一性原理研究中得出的数值符合得较好. 给出了两相Ce中如点缺陷形成能、表面能、层错能以及孪晶能等晶体缺陷形成能. 通过分析两相Ce的声子谱, 得出了不同温度下两相的晶格振动熵差, 其中在室温条件下约为0.67kB/atom. 还利用分子动力学模拟得出了该相变的等温线, 并且利用径向分布函数分析了相变前后两相的晶体结构, 确认了该相变为面心立方同构相变, 即Ce的α-γ相变. 由此表明, 本文的嵌入原子法势, 不仅可以分别合理地描述γ-Ce和α-Ce, 还可以反映γ-Ce和α-Ce两相之间的相变.
    Ce is a rare earth element in the periodic table. In the range of low temperature and low pressure, there are two face-centered-cubic (FCC) phases (α-Ce and γ-Ce) and a double-hexagonal-close-packed phase (β-Ce) for metallic Ce. At ambient temperature and about 0.7 GPa pressure, Ce undergoes γα phase transition with a volume shrink of 14%–17% discontinuously. In this paper, an embedded-atom method (EAM) potential compatible for α-Ce and γ-Ce was developed. This EAM potential has been employed to study several basic properties of cerium in these two FCC phases, such as equilibrium lattice constants, cohesive energies, and elastic constants. These results showed good accordance with experiments and first principle calculations. The lattice defects have been studied with the formation energy calculations of vacancies, interstitials, surfaces, stacking faults, and twinning defects in α-Ce and γ-Ce lattice. The lattice dynamics of α-Ce and γ-Ce have been analyzed using our EAM potential. The lattice vibrational entropy was calculated and plotted as functions of temperature for each phases. The vibrational entropy change across the α-γ phase transition showed to be ~0.67 kB per atom at ambient temperature. Using molecular dynamics simulation with our EAM potential, several isotherms and radial distribution functions were calculated. These isotherms and radial distribution functions demonstrate a first order phase transition between two FCC structures, corresponding to α-Ce and γ-Ce, with a critical point sets at Tc≈550 K and Pc≈1.21 GPa. Thus the newly developed EAM potential could provide a reasonable description of FCC Ce and its α-γ phase transition within the scale of classical molecular dynamics simulation.
      通信作者: 胡晓棉, hu_xiaomian@iapcm.ac.cn
      Corresponding author: Hu Xiao-Mian, hu_xiaomian@iapcm.ac.cn
    [1]

    郑海冰 2017 硕士学位论文 (哈尔滨: 哈尔滨工业大学

    Zheng H B 2017 M. S. Thesis (Harbin: Harbin Institute of Technoloty) (in Chinese)

    [2]

    金铭 2018 硕士学位论文 (哈尔滨: 哈尔滨工业大学

    Jin M 2018 M. S. Thesis (Harbin: Harbin Institute of Technoloty) (in Chinese)

    [3]

    林河成 2005 中国有色冶金 3 31Google Scholar

    Lin H C 2005 China Nonferrous Metallurgy 3 31Google Scholar

    [4]

    Koskenmaki D C, Gschneidner K A 1978 Handbook on the Physics and Chemistry of Rare Earths (Vol. 1) (Amsterdam: Elsevier North-Holland) pp337–377

    [5]

    潘昊, 胡晓棉, 吴子辉, 戴诚达, 吴强 2012 物理学报 61 206401Google Scholar

    Pan H, Hu X M, Wu Z H, Dai C D, Wu Q 2012 Acta Phys. Sin. 61 206401Google Scholar

    [6]

    Wang Y, Jr Hector L G, Zhang H, Shang S L, Chen L Q, Liu Z K 2008 Phys. Rev. B 78 104113Google Scholar

    [7]

    Decremps F, Belhadi L, Farber D L, Moore K T, Occelli F, Gauthier M, Polian A, Antonangeli D, Aracne-Ruddle C M, Amadon B 2011 Phys. Rev. Lett. 106 065701

    [8]

    Lipp M J, Jackson D, Cynn H, Aracne C, Evans W J, McMahan A K 2008 Phys. Rev. Lett. 101 165703Google Scholar

    [9]

    Johansson B 1974 Philos. Mag. 30 469Google Scholar

    [10]

    Allen J W, Martin R M 1982 Phys. Rev. Lett. 49 1106Google Scholar

    [11]

    Casadei M, Ren X, Rinke P, Rubio A, Scheffler M 2016 Phys. Rev. B 93 075153Google Scholar

    [12]

    Amadon B, Biermann S, Georges A, Aryasetiawan F 2006 Phys. Rev. Lett. 96 066402Google Scholar

    [13]

    Jeong I K, Darling T W, Graf M J, Proffen T, Heffner R H, Lee Y, Vogt T, Jorgensen J D 2004 Phys. Rev. Lett. 92 105702Google Scholar

    [14]

    El'kin V M, Kozlov E A, Kakshina E V, Moreva Y S 2006 Phys. Met. Metall. 101 208Google Scholar

    [15]

    El'kin V M, Mikhaylov V N, Petrovtsev A V, Cherne F J 2011 Phys. Rev. B 84 094120Google Scholar

    [16]

    Yelkin V M, Kozlov E A, Kakshina E V, Moreva Y S 2006 Shock Compression of Condensed Matter2005 Baltimore, Maryland 31 July–5 August, 2005 pp77–80

    [17]

    Casadei M, Ren X, Rinke P, Rubio A, Scheffler M 2012 Phys. Rev. Lett. 109 146402Google Scholar

    [18]

    Huang L, Chen CA 2007 J. Phys.Condens. Matter 19 476206Google Scholar

    [19]

    Krisch M, Farber D L, Xu R, Antonangeli D, Aracne C M, Beraud A, Chiang T-C, Zarestky J, Kim D Y, Isaev E I, Ahuja R, Johansson B 2011 Proc. Natl. Acad. Sci. U.S.A. 108 9342Google Scholar

    [20]

    Singh N, Singh S P 1990 Phys. Rev. B 42 1652Google Scholar

    [21]

    Hachiya K, Ito Y 1999 J. Phys. Condens. Matter 11 6543Google Scholar

    [22]

    Sheng H W, Kramer M J, Cadien A, Fujita T, Chen M W 2011 Phys. Rev. B 83 134118Google Scholar

    [23]

    Fu J, Zhao J 2013 Modell. Simul. Mater. Sci. Eng. 21 065003Google Scholar

    [24]

    Voter A F, Chen S P 1986 MRS Proceedings 82 175Google Scholar

    [25]

    Dupont V, Chen S P, Germann T C 2010 EPJ Web of Conferences Paris, France, May 24–28, 2010 p00009

    [26]

    Dupont V, Germann T C 2010 APS March Meeting Portland, Oregon, March 15–19, 2010 W30.00007

    [27]

    Rose J H, Smith J R, Guinea F, Ferrante J 1984 Phys. Rev. B 29 2963

    [28]

    Duff A I, Finnis M W, Maugis P, Thijsse B J, Sluiter M H F 2015 Comput. Phys. Commun. 196 439Google Scholar

    [29]

    Plimpton S 1995 J. Comput. Phys. 117 1Google Scholar

    [30]

    Eriksson O, Brooks M S S, Johansson B 1990 Phys. Rev. B 41 7311Google Scholar

    [31]

    Greiner J D, McMasters O D, Smith J F 1980 Scr. Metall. 14 989Google Scholar

    [32]

    Staun Olsen J, Gerward L, Dancausse J P, Gering E 1993 Physica B 190 92Google Scholar

    [33]

    Voronov F F, Goncharova V A, Stal'gorova O V 1979 Soviet Phys. JETP 49 687

    [34]

    Söderlind P, Eriksson O, Wills J M, Boring A M 1993 Phys. Rev. B 48 9306Google Scholar

    [35]

    Korzhavyi P A, Abrikosov I A, Johansson B, Ruban A V, Skriver H L 1999 Phys. Rev. B 59 11693Google Scholar

    [36]

    Press W H, Teukolsky S A, Vetterling W T, Flannery B P 2007 Numerical Recipes: The Art of Scientific Computing (3rd Ed.) (New York: Cambridge University Press) pp487–562

    [37]

    Östlin A, Di Marco I, Locht I L M, Lashley J C, Vitos L 2016 Phys. Rev. B 93 094103Google Scholar

    [38]

    van Swygenhoven H, Derlet P M, Frøseth A G 2004 Nat. Mater. 3 399Google Scholar

    [39]

    Hai S, Tadmor E B 2003 Acta Mater. 51 117Google Scholar

    [40]

    Tadmor E, Hai S 2003 J. Mech. Phys. Solids 51 765Google Scholar

    [41]

    Stassis C, Gould T, McMasters O D, Gschneidner K A, Nicklow R M 1979 Phys. Rev. B 19 5746Google Scholar

  • 图 1  Ce的EAM势 (a)对势函数和电子密度分布函数; (b)嵌入能函数

    Fig. 1.  EAM potential for Ce: (a) The pair function Φ(rij) and the density function f(rij); (b) the embedded function F(ρ).

    图 2  由本文EAM势得出的面心立方Ce的冷能曲线

    Fig. 2.  The cold energy of fcc Ce calculated with the newly developed EAM potential.

    图 3  面心立方Ce中的广义层错能和广义孪晶能 (a) α-Ce; (b) γ-Ce

    Fig. 3.  Generalized stacking fault energy and generalized twining fault energy curve for (a) α-Ce and (b) γ-Ce.

    图 4  面心立方Ce的声子态密度

    Fig. 4.  Phonon density of states for FCC Ce.

    图 5  α-Ce和γ-Ce两相的晶格振动熵Svib以及熵差ΔSvib (右下插图)随温度T的变化

    Fig. 5.  The vibrational entropy Svib and its change ΔSvib between two phases (the inset) plotted as functions of temperature.

    图 6  面心立方Ce的声子色散谱

    Fig. 6.  Phonon dispersion relations for FCC Ce.

    图 7  面心立方Ce的等温线 (a) NPT增压; (b) NPT减压; (c) NVT

    Fig. 7.  Isotherms of FCC Ce: (a) NPT, pressure increased; (b) NPT, pressure decreased; (c) NVT.

    图 8  Ce的γααγ相变路径

    Fig. 8.  The path of Ce γα and αγ phase transition.

    图 9  零压以及图7A, B两点状态的径向分布函数

    Fig. 9.  Radial distribution function for zero pressure (black), the A state (red), and the B state (blue) pointed in Fig.7.

    表 1  面心立方Ce的基本性质的EAM计算值与实验和第一性原理结果的比较

    Table 1.  EAM predicted properties of Ce lattice in comparison with experimental and ab initiodata.

    γ-Ceα-Ce
    实验第一性原理本文EAM实验第一性原理本文EAM
    a05.16[4]5.22[11]5.144.84[4] 4.90[30]4.63[11]4.81
    Ecoh/eV4.32[4]4.35[11]4.324.3[11]3.76[11]4.3255
    体弹模量/GPa18.18[31]28.3[11]16.7835.0[32], 16.94[33]37.0[34]37.00
    c11 /GPa26.01[31]23.0652.9[34]59.77
    c12 /GPa14.26[31]13.6429.1[34]25.62
    c44 /GPa17.30[31]17.6444.6[34]49.98
    剪切模量/GPa12.73[31]12.4717.26[33]36.82
    下载: 导出CSV

    表 2  γ-Ce and α-Ce中晶体缺陷的形成能

    Table 2.  Calculated formation energy of lattice defectsin γ-Ce and α-Ce.

    γ-Ceα-Ce
    之前的结果本文EAM之前的结果本文EAM
    Eif/eV3.3[22]1.932.97
    Evf/eV0.75[22], 2.02[23]0.851.15
    γ(100)/mJ·m–2697[22], 2140[23]391308
    γ(110)/mJ·m–2797[22], 2220[23]442390
    γ(111)/mJ·m–2586[22], 2190[23]297195
    γssf/mJ·m–2486[22],58[37], 16[37], –0.2[37]457301[37], 311[37], 369[37]734
    γusf/mJ·m–2501[22]543822
    γutf/mJ·m–212[22]7681167
    下载: 导出CSV
  • [1]

    郑海冰 2017 硕士学位论文 (哈尔滨: 哈尔滨工业大学

    Zheng H B 2017 M. S. Thesis (Harbin: Harbin Institute of Technoloty) (in Chinese)

    [2]

    金铭 2018 硕士学位论文 (哈尔滨: 哈尔滨工业大学

    Jin M 2018 M. S. Thesis (Harbin: Harbin Institute of Technoloty) (in Chinese)

    [3]

    林河成 2005 中国有色冶金 3 31Google Scholar

    Lin H C 2005 China Nonferrous Metallurgy 3 31Google Scholar

    [4]

    Koskenmaki D C, Gschneidner K A 1978 Handbook on the Physics and Chemistry of Rare Earths (Vol. 1) (Amsterdam: Elsevier North-Holland) pp337–377

    [5]

    潘昊, 胡晓棉, 吴子辉, 戴诚达, 吴强 2012 物理学报 61 206401Google Scholar

    Pan H, Hu X M, Wu Z H, Dai C D, Wu Q 2012 Acta Phys. Sin. 61 206401Google Scholar

    [6]

    Wang Y, Jr Hector L G, Zhang H, Shang S L, Chen L Q, Liu Z K 2008 Phys. Rev. B 78 104113Google Scholar

    [7]

    Decremps F, Belhadi L, Farber D L, Moore K T, Occelli F, Gauthier M, Polian A, Antonangeli D, Aracne-Ruddle C M, Amadon B 2011 Phys. Rev. Lett. 106 065701

    [8]

    Lipp M J, Jackson D, Cynn H, Aracne C, Evans W J, McMahan A K 2008 Phys. Rev. Lett. 101 165703Google Scholar

    [9]

    Johansson B 1974 Philos. Mag. 30 469Google Scholar

    [10]

    Allen J W, Martin R M 1982 Phys. Rev. Lett. 49 1106Google Scholar

    [11]

    Casadei M, Ren X, Rinke P, Rubio A, Scheffler M 2016 Phys. Rev. B 93 075153Google Scholar

    [12]

    Amadon B, Biermann S, Georges A, Aryasetiawan F 2006 Phys. Rev. Lett. 96 066402Google Scholar

    [13]

    Jeong I K, Darling T W, Graf M J, Proffen T, Heffner R H, Lee Y, Vogt T, Jorgensen J D 2004 Phys. Rev. Lett. 92 105702Google Scholar

    [14]

    El'kin V M, Kozlov E A, Kakshina E V, Moreva Y S 2006 Phys. Met. Metall. 101 208Google Scholar

    [15]

    El'kin V M, Mikhaylov V N, Petrovtsev A V, Cherne F J 2011 Phys. Rev. B 84 094120Google Scholar

    [16]

    Yelkin V M, Kozlov E A, Kakshina E V, Moreva Y S 2006 Shock Compression of Condensed Matter2005 Baltimore, Maryland 31 July–5 August, 2005 pp77–80

    [17]

    Casadei M, Ren X, Rinke P, Rubio A, Scheffler M 2012 Phys. Rev. Lett. 109 146402Google Scholar

    [18]

    Huang L, Chen CA 2007 J. Phys.Condens. Matter 19 476206Google Scholar

    [19]

    Krisch M, Farber D L, Xu R, Antonangeli D, Aracne C M, Beraud A, Chiang T-C, Zarestky J, Kim D Y, Isaev E I, Ahuja R, Johansson B 2011 Proc. Natl. Acad. Sci. U.S.A. 108 9342Google Scholar

    [20]

    Singh N, Singh S P 1990 Phys. Rev. B 42 1652Google Scholar

    [21]

    Hachiya K, Ito Y 1999 J. Phys. Condens. Matter 11 6543Google Scholar

    [22]

    Sheng H W, Kramer M J, Cadien A, Fujita T, Chen M W 2011 Phys. Rev. B 83 134118Google Scholar

    [23]

    Fu J, Zhao J 2013 Modell. Simul. Mater. Sci. Eng. 21 065003Google Scholar

    [24]

    Voter A F, Chen S P 1986 MRS Proceedings 82 175Google Scholar

    [25]

    Dupont V, Chen S P, Germann T C 2010 EPJ Web of Conferences Paris, France, May 24–28, 2010 p00009

    [26]

    Dupont V, Germann T C 2010 APS March Meeting Portland, Oregon, March 15–19, 2010 W30.00007

    [27]

    Rose J H, Smith J R, Guinea F, Ferrante J 1984 Phys. Rev. B 29 2963

    [28]

    Duff A I, Finnis M W, Maugis P, Thijsse B J, Sluiter M H F 2015 Comput. Phys. Commun. 196 439Google Scholar

    [29]

    Plimpton S 1995 J. Comput. Phys. 117 1Google Scholar

    [30]

    Eriksson O, Brooks M S S, Johansson B 1990 Phys. Rev. B 41 7311Google Scholar

    [31]

    Greiner J D, McMasters O D, Smith J F 1980 Scr. Metall. 14 989Google Scholar

    [32]

    Staun Olsen J, Gerward L, Dancausse J P, Gering E 1993 Physica B 190 92Google Scholar

    [33]

    Voronov F F, Goncharova V A, Stal'gorova O V 1979 Soviet Phys. JETP 49 687

    [34]

    Söderlind P, Eriksson O, Wills J M, Boring A M 1993 Phys. Rev. B 48 9306Google Scholar

    [35]

    Korzhavyi P A, Abrikosov I A, Johansson B, Ruban A V, Skriver H L 1999 Phys. Rev. B 59 11693Google Scholar

    [36]

    Press W H, Teukolsky S A, Vetterling W T, Flannery B P 2007 Numerical Recipes: The Art of Scientific Computing (3rd Ed.) (New York: Cambridge University Press) pp487–562

    [37]

    Östlin A, Di Marco I, Locht I L M, Lashley J C, Vitos L 2016 Phys. Rev. B 93 094103Google Scholar

    [38]

    van Swygenhoven H, Derlet P M, Frøseth A G 2004 Nat. Mater. 3 399Google Scholar

    [39]

    Hai S, Tadmor E B 2003 Acta Mater. 51 117Google Scholar

    [40]

    Tadmor E, Hai S 2003 J. Mech. Phys. Solids 51 765Google Scholar

    [41]

    Stassis C, Gould T, McMasters O D, Gschneidner K A, Nicklow R M 1979 Phys. Rev. B 19 5746Google Scholar

  • [1] 张晓军, 王安祥, 严祥安, 陈长乐. 改进分析型嵌入原子法在W(100)表面声子谱中的应用. 物理学报, 2020, 69(7): 076301. doi: 10.7498/aps.69.20191910
    [2] 第伍旻杰, 胡晓棉. 单晶Ce冲击相变的分子动力学模拟. 物理学报, 2020, 69(11): 116202. doi: 10.7498/aps.69.20200323
    [3] 田晓林, 赵宇宏, 田晋忠, 侯华. 原子间相互作用势对中Al浓度Ni75AlxV25-x合金沉淀序列的影响. 物理学报, 2018, 67(23): 230201. doi: 10.7498/aps.67.20181366
    [4] 孙素蓉, 王海兴. 惰性气体原子间相互作用势比较研究. 物理学报, 2015, 64(14): 143401. doi: 10.7498/aps.64.143401
    [5] 赵健东, 辛洁. 高激发态原子间的范德瓦尔斯相互作用. 物理学报, 2014, 63(13): 133201. doi: 10.7498/aps.63.133201
    [6] 王晓璐, 令狐荣锋, 宋晓书, 吕兵, 杨向东. 氦原子与卤族氢化物分子相互作用势的研究. 物理学报, 2013, 62(16): 163101. doi: 10.7498/aps.62.163101
    [7] 令狐荣锋, 徐梅, 吕兵, 宋晓书, 杨向东. He原子与N2分子相互作用势的理论研究. 物理学报, 2013, 62(1): 013103. doi: 10.7498/aps.62.013103
    [8] 安荣, 刘威, 王春青, 田艳红. Sn的分析型键级势. 物理学报, 2013, 62(12): 128101. doi: 10.7498/aps.62.128101
    [9] 潘昊, 胡晓棉, 吴子辉, 戴诚达, 吴强. 铈低压冲击相变数值模拟研究. 物理学报, 2012, 61(20): 206401. doi: 10.7498/aps.61.206401
    [10] 敖冰云, 汪小琳, 陈丕恒, 史鹏, 胡望宇, 杨剑瑜. 嵌入原子法计算金属钚中点缺陷的能量. 物理学报, 2010, 59(7): 4818-4825. doi: 10.7498/aps.59.4818
    [11] 陈怡, 申江. NaZn13型Fe基化合物的结构和热力学性质研究. 物理学报, 2009, 58(13): 141-S145. doi: 10.7498/aps.58.141
    [12] 陈怡, 申江. 稀土金属间化合物RFe2Zn20-xInx的结构属性研究. 物理学报, 2009, 58(13): 146-S150. doi: 10.7498/aps.58.146
    [13] 江慧丰, 张青川, 陈学东, 范志超, 陈忠家, 伍小平. 位错与溶质原子间动态相互作用的数值模拟研究. 物理学报, 2007, 56(6): 3388-3392. doi: 10.7498/aps.56.3388
    [14] 王永亮, 张 超, 唐 鑫, 张庆瑜. 表面Cu原子间相互作用对Cu(001)表面跳跃扩散行为的影响. 物理学报, 2006, 55(8): 4214-4220. doi: 10.7498/aps.55.4214
    [15] 李 明, 孙久勋. 原子间相互作用对光场和原子激光压缩性质的影响. 物理学报, 2006, 55(6): 2702-2707. doi: 10.7498/aps.55.2702
    [16] 周明, 黄春佳. 原子间相互作用对原子激光压缩性质的影响. 物理学报, 2004, 53(1): 54-57. doi: 10.7498/aps.53.54
    [17] 张建民, 徐可为, 马 飞. 用改进嵌入原子法计算Cu晶体的表面能. 物理学报, 2003, 52(8): 1993-1999. doi: 10.7498/aps.52.1993
    [18] 李 泌. 铁的原子间相互作用及声子谱. 物理学报, 2000, 49(9): 1692-1695. doi: 10.7498/aps.49.1692
    [19] 万钧, 申三国, 范希庆. Cu表面弛豫和自扩散机制的修正嵌入原子法模拟. 物理学报, 1997, 46(6): 1161-1167. doi: 10.7498/aps.46.1161
    [20] 刘家瑞, 张建华, 李毅. 晶格原子间非简谐相互作用对离子注入的影响. 物理学报, 1989, 38(9): 1400-1405. doi: 10.7498/aps.38.1400
计量
  • 文章访问数:  5600
  • PDF下载量:  81
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-06-06
  • 修回日期:  2019-08-09
  • 上网日期:  2019-10-01
  • 刊出日期:  2019-10-20

面心立方Ce同构相变的分子动力学模拟

  • 1. 中国工程物理研究院研究生院, 北京 100088
  • 2. 北京应用物理与计算数学研究所, 计算物理国家重点实验室, 北京 100088
  • 通信作者: 胡晓棉, hu_xiaomian@iapcm.ac.cn

摘要: 基于嵌入原子法, 本文给出了一个金属Ce原子间的相互作用势. 利用该势分别计算了γ-Ce和α-Ce的晶格常数、结合能、弹性常数, 计算结果与实验或第一性原理研究中得出的数值符合得较好. 给出了两相Ce中如点缺陷形成能、表面能、层错能以及孪晶能等晶体缺陷形成能. 通过分析两相Ce的声子谱, 得出了不同温度下两相的晶格振动熵差, 其中在室温条件下约为0.67kB/atom. 还利用分子动力学模拟得出了该相变的等温线, 并且利用径向分布函数分析了相变前后两相的晶体结构, 确认了该相变为面心立方同构相变, 即Ce的α-γ相变. 由此表明, 本文的嵌入原子法势, 不仅可以分别合理地描述γ-Ce和α-Ce, 还可以反映γ-Ce和α-Ce两相之间的相变.

English Abstract

    • 铈(cerium, Ce)是一种稀土金属元素, 其f电子处于局域与非局域的转变点上, 具有低熔点、非对称晶体结构、多种同素异形体共存、发生相变伴随体积变化等物理、化学性质[1,2]. Ce及其稀土合金已广泛应用于钢铁、有色金属及其合金、发火合金、稀土永磁合金、储氢材料等诸多领域[3].

      在较低的温度和压强条件下, 金属Ce有两种面心立方(FCC)相(α相和γ相)以及一种双密排六方相(β相)[4]. 在室温和0.7GPa压缩条件下Ce会发生γα同构相变[4-6], 晶格常数从5.16 Å突变到4.84 Å, 同时体积减小14%—17%. α-γ相变的临界点(C.P.)压强PC = 1.44—1.96 GPa[4,7,8], 温度TC = 460—600 K[4,7,8].

      作为Ce的α-γ相变的机制, 从电子结构来看, Mott相变和Kondo体积坍缩(KVC)两种理论模型长期争论不下[4,9,10]. 也有研究认为, 这两种模型分别反映了Ce的α-γ相变的不同侧面[11]. 除此之外, Amadon等[12]认为在室温条件下熵的效应不可忽略, Ce的α-γ相变是由熵驱动的. 而熵的贡献可分为电子的和晶格振动的两个自由度. 利用高压高分辨率同步X射线粉末衍射实验, Jeong等[13]测得Ce的α-γ相变过程中, 晶格振动的熵变(ΔSvibα-γ ≈ (0.75 ± 0.15)kB/atom)约占总熵变(约1.54kB/atom)的一半. Wang等[6]利用第一性原理计算得出的总熵变约为2.86kB/atom, 其中晶格振动的熵变约为0.94kB/atom. 这些都表明晶格振动的熵变在总熵变中占据相当可观的一部分.

      对于Ce的α-γ相变, 以往的模拟计算工作中主要分宏观数值模拟和第一性原理电子结构计算两方面. 基于多相物态方程的宏观数值模拟[14-16]无法给出相变过程中微尺度结构的细节; 而相关的第一性原理计算[6,11,17-19]受制于其时间和空间尺度, 只能对此相变进行均匀、平衡态的分析. 而非平衡态分子动力学模拟因其特点, 可以在微、纳米尺度提供更多的信息, 如相变过程中微结构的演化, 从而可以将离散原子尺度与宏观尺度沟通起来.

      在分子动力学模拟中, 分子间相互作用势必不可少. 已有工作中针对FCC结构Ce的经验势不多: 例如Singh和Singh[20]将两体作用的杂化近自由电子紧束缚成键模型用在一些FCC结构的稀土金属(包括Ce), 得出的结果在声子谱和弹性常数等方面与实验符合较好; Hachiya和Ito[21]在此基础上利用近自由电子键角相关紧束缚成键模型(最早用于过渡金属及其合金)进行了改进; Sheng等[22]利用势能面拟合法, 得出了包括Ce的14中FCC结构金属的嵌入原子法(EAM)势; Fu和Zhao[23]拟合了FCC镧(La)和Ce的Gupta势. 虽然这些经验作用势都与Ce的实验和第一性原理的数据符合得很好, 但只涉及了γ-Ce, 而与α-Ce无关, 因此无法用于Ce的α-γ相变的分子动力学模拟. Voter和Chen[24]通过对Voter-Chen EAM进行修正, 提高两相之间的势垒(0.142 eV/atom)来维持各相的稳定性, 提出过一种能够描述Ce的α-γ相变的经验势[25,26], 不过该势的具体描述至今未见发表. 为了能够对Ce的α-γ相变进行大规模分子动力学模拟,本文展示了一种新的FCC结构Ce的EAM势, 该势可对γ-Ce, α-Ce以及Ce的α-γ相变进行合理描述.

    • α-Ce和γ-Ce的冷能曲线可以分别用Rose物态方程来表示[25,27], 也就是说FCC结构Ce的冷能曲线应当是一个分段函数: 晶格常数较大的部分为γ-Ce的Rose曲线, 晶格常数较小的部分则为α-Ce的Rose曲线. 两段各有一个极小值点对应于各相平衡位置的晶格常数, 两段曲线交叉于两极小值点中间形成一个势垒.

      在嵌入原子法中, 单类型原子系统的总势能Etot可以表示为

      $ \begin{split} & {E_{{\rm{tot}}}} = \frac{1}{2}\mathop \sum \limits_{i,j} \varPhi \left( {{r_{ij}}} \right) + \mathop \sum \limits_i F\left( {{\rho _i}} \right) ,\\ &{\rho _i} = \mathop \sum \limits_{j \ne i} f\left( {{r_{ij}}} \right),\end{split}$

      其中Φ(rij)为原子ij之间的对势, F为原子i处于电子云密度ρi的嵌入能, f(rij)为原子j在原子i处的电子密度分布函数.

      在文献[28]中, 嵌入能F和电子密度ρ的关系表示为

      $ F\left( \rho \right) = a{\rho ^{1/2}} + b{\rho ^2} + c{\rho ^3}, $

      式中a, b, c为待定参数. Φ(rij)和f(rij)的解析形式为立方项求和:

      $ \begin{split} & \varPhi \left( r \right) = \mathop \sum \limits_n {a_n}{\left( {{r_n} - r} \right)^3}\varTheta \left( {{r_n} - r} \right),\\ &f\left( r \right) = \mathop \sum \limits_n {b_n}{\left( {{s_n} - r} \right)^3}\varTheta \left( {{s_n} - r} \right). \end{split} $

      可以根据精度要求选取n的上限, an, rn, bn, sn为待定参数; Θ(x)为Heaviside函数, 若x ≤ 0, Θ(x) = 0; 若x > 0, Θ(x) = 1. 已经有研究将其用于拟合面心立方金属的原子间相互作用势. 在本文工作中先对Φ(rij)和f(rij)分别只取一项, 即

      $ \varPhi \left( r \right) = {a^\gamma }{\left( {{r^\gamma } - r} \right)^3}\varTheta \left( {{r^\gamma } - r} \right), $

      $ f\left( r \right) = {b^\gamma }{\left( {{s^\gamma } - r} \right)^3}\varTheta \left( {{s^\gamma } - r} \right), $

      且令bγ = 1 Å–3, sγ = 6.6 Å (此即相互作用势的截断距离). 利用最小二乘法, 以γ-Ce的Rose物态方程曲线和弹性常数为目标对(2)式—(4)式中的参数a, b, c, aγ, rγ进行拟合. 由此得出的γ-Ce的冷能曲线与α-Ce的Rose物态方程曲线相交于r = 3.498 Å, 于是对Φ(r), f(r)和F(ρ)函数应当采用分段函数的形式(分段点为r = 3.498 Å, ρ = 389.094). 由于势函数的长程部分(r > 3.498 Å) 已经确定, 接下来只能调整短程的部分(r < 3.498 Å). 在r < 3.498 Å的范围内仍采用相似的形式, 同时为了保证势函数的连续性, 则有

      $ \varPhi \left( r \right) = {a^\alpha }{\left( {{r^\alpha } - r} \right)^3}\varTheta \left( {{r^\alpha } - r} \right) + {\varPhi _0}, $

      $ f\left( r \right) = {b^\alpha }{\left( {{s^\alpha } - r} \right)^3}\varTheta \left( {{s^\alpha } - r} \right) + {f_0}. $

      嵌入函数F(ρ)则通过嵌入能F (即总势能与对势能部分之差)与电子云密度ρ的对应关系利用立方样条插值得出. (5)式和(6)式中的参数aα, rα, Φ0, bα, sα, f0参照α-Ce的Rose物态方程曲线和弹性常数, 并考虑α-γ相变的温度和压强条件来确定. 最终的对势函数Φ(r)、电子密度函数f(r)、嵌入能函数F(ρ)如图1(a)图1(b)所示.

      图  1  Ce的EAM势 (a)对势函数和电子密度分布函数; (b)嵌入能函数

      Figure 1.  EAM potential for Ce: (a) The pair function Φ(rij) and the density function f(rij); (b) the embedded function F(ρ).

    • 为了验证第2节中所得的EAM势, 本文采用了LAMMPS[29]对面心立方Ce晶体的一些性质以及α-γ相变过程进行了分子动力学模拟计算.

    • 图2展示了由EAM势计算得到的Ce的冷能曲线. 该曲线分为两段, 每段各有一个极小点, 分别对应于α-Ce (a = 4.81 Å)和γ-Ce (a = 5.14 Å). 可以看出, 在两个极小值点中间存在一个高度为15.3—20.8 meV/atom的势垒, 远小于Chen等[25]的经验势中的对应值(0.142 eV/atom). 表1中展示了由EAM计算得出的α-Ce和γ-Ce的一些基本性质, 包括晶格常数、结合能、弹性常数等, 与实验或第一性原理计算的结果符合得很好.

      图  2  由本文EAM势得出的面心立方Ce的冷能曲线

      Figure 2.  The cold energy of fcc Ce calculated with the newly developed EAM potential.

      γ-Ceα-Ce
      实验第一性原理本文EAM实验第一性原理本文EAM
      a05.16[4]5.22[11]5.144.84[4] 4.90[30]4.63[11]4.81
      Ecoh/eV4.32[4]4.35[11]4.324.3[11]3.76[11]4.3255
      体弹模量/GPa18.18[31]28.3[11]16.7835.0[32], 16.94[33]37.0[34]37.00
      c11 /GPa26.01[31]23.0652.9[34]59.77
      c12 /GPa14.26[31]13.6429.1[34]25.62
      c44 /GPa17.30[31]17.6444.6[34]49.98
      剪切模量/GPa12.73[31]12.4717.26[33]36.82

      表 1  面心立方Ce的基本性质的EAM计算值与实验和第一性原理结果的比较

      Table 1.  EAM predicted properties of Ce lattice in comparison with experimental and ab initiodata.

    • 为了定量地研究在Ce中的点缺陷, 采用本文的EAM势在10 × 10 × 10 的FCC超胞中计算了α-Ce和γ-Ce中的空位和填隙原子的形成能. 晶体点缺陷的形成能(Ef)的定义为

      $ {E_{\rm{f}}} = {E_{{\rm{def}}}} - \left( {N \pm 1} \right){E_{{\rm{coh}}}}, $

      其中Edef为包含单个点缺陷超胞的最小化总能量, N为无缺陷超胞的原子数(在此N = 4000). 于是正号对应于填隙原子, 负号对应于空位. Ecoh为无缺陷晶体中单个原子的结合能.

      对于γ-Ce, 填隙原子和空位的形成能分别为1.93和0.85 eV, 分别对应于无缺陷单晶中单个原子结合能的44.7%和19.7%; α-Ce中的填隙原子和空位的形成能分别为2.97和1.15 eV, 分别对应于无缺陷单晶中单个原子结合能的68.7%和26.6%. 以往的研究中所得的FCC过渡金属的空位能约为1—3 eV[35], 这里所得的填隙原子和空位的形成能处在与之可比的范围内, 可以认为是合理的.

    • 采用本文的EAM势同时通过Polak-Ribiere共轭梯度法[36]对表面构型进行充分弛豫, 以此计算了α-Ce和γ-Ce的三种低密勒指数(即(100), (110)和(111))晶面的表面形成能, 计算结果列在了表2中. 对于γ-Ce本文中计算所得的表面能明显低于以往经验势的计算结果[22,23]. 不过两相Ce中的表面能的大小顺序都符合γ(111) < γ(100) < γ(110), 表明其中(111)表面最稳定, 这方面与Sheng等[22]以及Fu和Zhao[23]γ-Ce的计算相结果相一致. 比较α-Ce和γ-Ce相同晶向的表面能, 则有γα < γγ, 如对于(100)晶面γα(100) < γγ(100), 另外对于γ(110)γ(111)亦同. 到目前尚未见关于金属Ce表面能实验测量结果的报道, 本文的表面能数值可作为未来相关测量实验的参考.

      γ-Ceα-Ce
      之前的结果本文EAM之前的结果本文EAM
      Eif/eV3.3[22]1.932.97
      Evf/eV0.75[22], 2.02[23]0.851.15
      γ(100)/mJ·m–2697[22], 2140[23]391308
      γ(110)/mJ·m–2797[22], 2220[23]442390
      γ(111)/mJ·m–2586[22], 2190[23]297195
      γssf/mJ·m–2486[22],58[37], 16[37], –0.2[37]457301[37], 311[37], 369[37]734
      γusf/mJ·m–2501[22]543822
      γutf/mJ·m–212[22]7681167

      表 2  γ-Ce and α-Ce中晶体缺陷的形成能

      Table 2.  Calculated formation energy of lattice defectsin γ-Ce and α-Ce.

    • 材料的稳定和非稳定层错能被认为是决定其力学行为(如金属的塑性变形)的重要物理量. 利用本文的EAM势计算得出了α-Ce和γ-Ce中堆垛层错和(111)孪晶缺陷的广义面缺陷能曲线. 广义层错能曲线可通过对无缺陷单晶沿(111)晶面进行刚性剪切得到, 而广义孪晶缺陷能曲线则可通过对预置层错的无缺陷单晶进行刚性剪切得出. 最终α-Ce和γ-Ce中的广义层错能曲线和广义孪晶能曲线分别如图3(a)图3(b)所示.

      图  3  面心立方Ce中的广义层错能和广义孪晶能 (a) α-Ce; (b) γ-Ce

      Figure 3.  Generalized stacking fault energy and generalized twining fault energy curve for (a) α-Ce and (b) γ-Ce.

      由此可得γ-Ce的非稳定层错能γusf = 0.543 J·m–2, 稳定层错能γssf = 0.457 J·m–2 (与Sheng等[22]利用EAM计算的结果相近); 而在α-Ce中γusf = 0.821 J·m–2, γssf = 0.734 J·m–2. 虽然此处的结果与Östlin等[37]的第一性原理计算结果相差较大, 但是对比α-Ce和γ-Ce两相的稳定层错能和非稳定层错能均满足γγsf < γαsf. Swygenhoven等[38]认为要了解FCC晶体变形的机制仅有非稳定层错能或者稳定层错能都是不全面的, 还应当比较稳定层错能与非稳定层错能的比值γssf/γusf. 该比值接近于1的时候倾向于越过整个势垒产生全位错; 反之若该比值较小则倾向于产生偏位错(作为参考, 在Al中该比值为0.97, 在Ni中为0.55, 在Cu中为0.13). 在γ-Ce中γssf/γusf = 0.842, α-Ce中γssfusf = 0.875, 数值较高, 与Al中的比值比较接近. 因此可以推测本文的EAM势在MD模拟的FCC结构Ce塑性变形中有可能观察到全位错产生.

      除位错激活与滑移之外, 形变孪晶是FCC结构晶体的另一种形变机制. 由本文EAM势计算得出的非稳定孪晶能γutf, 在γ-Ce中γutf = 0.768 J·m–2, 在α-Ce中γutf = 1.167 J·m–2. Tadmor和Hai[39,40]发现形成形变孪晶的趋势与比值γutf/γusf相关, 该比值越高越不容易产生形变孪晶. 对于γ-Ce, γutf/γusf = 1.414; 对于α-Ce, γutf/γusf = 1.420, 均接近且高于其在Al中的比值(约1.3). 由此可推测利用本文EAM势进行MD模拟很难出现形变孪晶.

    • 图4所示为温度T = 300 K条件下晶格常数分别等于4.89 Å (α-Ce, 压强约0.7 GPa)和 5.16 Å (γ-Ce, 零压)状态的声子态密度. 两相的态密度有显著差异: 1)“肩部”由1.21 THz (γ-Ce) 移动到了 1.43 THz (α-Ce); 2)后续增长的尖峰由1.81 THz (γ-Ce) 移动到了 2.31 THz (α-Ce); 3)再之后的峰由2.58 THz (γ-Ce)移动到了3.36 THz (α-Ce).

      图  4  面心立方Ce的声子态密度

      Figure 4.  Phonon density of states for FCC Ce.

      Ce的α-γ相变被认为是由熵驱动的, 也就是说熵在其中扮演重要位置. 利用图4中的态密度, 可以分别计算出两种状态的晶格振动的熵(如图5). 图5中右下角虚线图为两种状态的晶格振动熵的差值, 温度T = 300 K时的熵差为ΔSvib = 0.67kB/atom, 此与Jeong等[13]的结果(0.75 ± 0.15)kB/atom相近. 然而由于经典分子动力学模拟的特性, 无法给出电子结构自由度对熵的贡献, 在此对Ce的α-γ相变过程的经典分子动力学模拟中只能将固体系统的熵全部归于晶格振动自由度. 于是本文工作中分子动力学模拟的Ce的α-γ相变过程的总熵变等于晶格振动熵变, 即ΔS = ΔSvib = 0.67kB/atom (此值较实验值1.54kB/atom少约0.87 kB/atom).

      图  5  α-Ce和γ-Ce两相的晶格振动熵Svib以及熵差ΔSvib (右下插图)随温度T的变化

      Figure 5.  The vibrational entropy Svib and its change ΔSvib between two phases (the inset) plotted as functions of temperature.

    • 利用声子色散曲线可以更严格地检验经验作用势. 图6所示为300 K温度条件下计算利用本文EAM计算所得α-Ce (a = 4.89 Å, 黑实线) 和γ-Ce (a = 5.16 Å, 蓝色点虚线) 的声子色散关系谱. 两相均有三支声学波, 其中一支为纵波, 另两支为横波. 之前的实验和理论研究中发现对于γ-Ce在 [ζζζ] 方向横波的L点处存在反常下沉(软模)[6,19,41]. 而本文计算中无论α-Ce还是γ-Ce均未找到此反常下沉.

      图  6  面心立方Ce的声子色散谱

      Figure 6.  Phonon dispersion relations for FCC Ce.

    • 利用本文中的EAM势并通过分子动力学模拟研究了Ce的α-γ相变. 对初始为FCC结构的Ce进行静水压(各向同性)加载, 计算所得的压强-体积等温线如图7(a)图7(c)所示, 其中图7(a)图7(b)的等温线由NPT系综分别增、减压强得出, 图7(c)的等温线由NVT系综弛豫得出. 如在温度为T = 300 K的条件下, 增大压强到0.80 GPa时发生相变, 体积突变减小2.68 Å3/atom; 反之, 减小压强到0.59 GPa时发生逆相变, 体积突变增大2.83 Å3/atom. 由此看出300 K温度条件下存在明显的(0.21 GPa)迟滞.当温度升高至约550 K时, 图7(a)图7(b)中的等温线变为光滑曲线, 体积随压强连续变化. 对应图7(c)中的极小值点与极大值点重合于拐点处, 此时拐点的斜率为零. 当达到更高的温度, 如600, 700和800 K时, 拐点的斜率变为负值. 由此得出临界点温度Tc ≈ 550 K, 压强Pc ≈ 1.21 GPa. 根据图7(a)图7(b)绘出的静水压加载下该相变的温度-压强相变路径如图8所示, 此可视为该相变的P-T相图.

      图  7  面心立方Ce的等温线 (a) NPT增压; (b) NPT减压; (c) NVT

      Figure 7.  Isotherms of FCC Ce: (a) NPT, pressure increased; (b) NPT, pressure decreased; (c) NVT.

      图  8  Ce的γααγ相变路径

      Figure 8.  The path of Ce γα and αγ phase transition.

      为了确认该相变否是为Ce的α-γ同构相变, 利用径向分布函数对相变前后的晶格结构进行了分析. 图9中比较了温度为300 K时不同压强加载下的径向分布函数: 零压加载下的前三个峰值对应的间距值分别为3.62, 5.16 和6.30 Å; 增压加载当压强增大至0.7 GPa (图7A点所指状态) 时对应的这三个间距值分别为3.56, 5.06 和6.18 Å; 减压加载当压强降至0.7 GPa时 (图7B点所指状态) 对应的这三个间距值分别为3.42, 4.89 和5.98 Å. 这三种不同加载情况下的径向分布函数均符合FCC结构的最近三阶近邻间距之间的关系. 于是可以认为相变前后的两相均为FCC结构, 以上模拟的相变正是Ce的α-γ同构相变, 其中晶格常数较大的一相对应于γ-Ce, 而晶格常数较小的一相对应于α-Ce. 也就是说, 本文的EAM势可以用于FCC结构Ce的α-γ同构相变的分子动力学模拟, 对α-γ同构相变进行多尺度的研究.

      图  9  零压以及图7A, B两点状态的径向分布函数

      Figure 9.  Radial distribution function for zero pressure (black), the A state (red), and the B state (blue) pointed in Fig.7.

    • 根据Rose物态方程还有弹性常数等实验以及第一性原理计算数据, 本文给出了FCC结构Ce的EAM势. 利用该EAM势对两种FCC结构Ce进行了定量的讨论. 计算得出的两相Ce的结合能、弹性常数等一些基本性质, 以及晶体缺陷和声子谱, 与实验和第一性原理计算的结果符合较好. 此外该EAM势可以用于FCC结构Ce的α-γ同构相变的分子动力学模拟, 并且通过MD模拟得出的P-T相图也与实验相符合.

      需要注意的是, 以往的研究中表明Ce的α-γ相变是熵驱动的(特别是当处在室温条件下), 其中包括了电子自由度和晶格振动自由度两个方面. 但由于分子动力学模拟的局限性, 无法体现电子自由度在其中的贡献, 只能给出晶格振动方面的影响.本文中计算得出的常温下两相的晶格振动熵差约为0.67kB/atom.

      本文的EAM势主要针对的是FCC结构的α-Ce和γ-Ce以及Ce的α-γ同构相变, 未考虑到其他的因素. 如同样处于较低的温度和压强条件下的β相, 以及温度较高时的δ相(体心立方结构)和液态相. 对此需要在后续的工作中考虑并改进.

参考文献 (41)

目录

    /

    返回文章
    返回