氦热力学性质的分子动力学研究

第23卷第6期强激光与粒束Vol 23, No 62011年6月HIGH POWER LASER AND PARTICLE BEAMSJun.,2011文章编号:1001-4322(2011)06-1649-04氦热力学性质的分子动力学研究何以广1:2,王钊2,梁晶2,高爽2,田宝贤2,汤秀章2(1,清华大学工程物理系,北京100084;2.中国原子能科学研究院高功率准分子激光实验室,北京102413)摘要:利用经典分子动力学和第一性原理分子动力学,研究了氦在高压下的熔化曲线、状态方程和非金属-金属转变。得到了氦在温度小于4.5eV、密度0.3~5.0g/cm3范围内的状态方程,并把氮的熔化曲线的压强范围拓展到了50GPa氦的能隙宽度曲线表明,温度大大降低了氦的金属化密度。关键词:状态方程;分子动力学;非金属-金属转变;熔化曲线中图分类号:O521文献标志码doi:10.3788/ HPLPB20112306.1649氮是宇宙中含量最丰富的元素之一,它和氢组成的混合物是构成宇宙中巨型气体行星的主要物质11。为了描述这些行星的结构和演化,需要了解氦的熔化曲线状态方程等热力学数据27。实验上,利用静态高压技术,可测量氦的熔化曲线,但压强范围受到静态高压技术的限制。理论上,在低压区可以采用经验势函数和经典分子动力学得到氦的熔化温度。但在高压区,经验势函数不再适用。氦的状态方程也得到了广泛的研究。最近的冲击波压缩实验表明,氦的单次冲击波压缩率可达到6.00。理论上,在低压区可以用自由能最小原理得到氮的状态方程1,但这种方法不适用高压范围。第一性原理分子动力学从第一性原理出发,不需要任何经验参数,可以突破实验上的压强范围,特别适合高压物性的研究。因此本文用第一性原理分子动力学研究了氮在高压下的热力学性质。1计算理论和方法1.1第一性原理分子动力学在第一性原理分子动力学中,电子作量子力学处理离子作经典处理。因此核心的问题是如何得到电子的波函数。在本文中,我们用密度泛函理论(DFT)得到了电子的波函数。密度泛函理论以 Hohenberg和Kohn的理论为基础213),对于N个电子的体系,方程如下[-A-vi+Vu(r)]w, (r)=e,,(r)(1)[p(r)](r)Zierr-R,+Vxc Lp(r)]式中:V←m是有效势能;v(r)是波函数;p(r)是电子密度;Vxc是交换相关势;r,R均为粒子位置;e为电子电荷;z为原子序数。在本文中需要考虑温度的影响,因此在第一性原理分子动力学过程中,我们采用了 Mermin有限温度密度泛函理论141)。在计算中,利用了投影子缀加平面波膺势(PAW)和广义梯度近似交换相关势(GGA)1·7,截断能量Ecm为900eV。采用了108个原子的超胞近似,并用了周期性边界条件。温度用Nos恒温器控制,let算法,时间步长为0.4~2.0fs,密度越大,步长越短。体系达到平衡后,体系继续运行1000~200r电子温度用电子能带的 Fermi分布确定。空间积分用布里渊的T点采样。在分子动力学过程中,利用了V1.2经典分子动力学在经典分子动力学中,氦原子间的相互作用势能采用exp6有效势函数φ(r)6exp[a(1-)]-a()6}3)由于→→0时,exp6势函数没有物理意义,因此当r≤W时,exp6势函数用φ=Aexp(-Br)代替。其中,a=中国煤化工收稿日期:2010-0510;修订日期:2011-0428CNMHG基金項目:国网家高技术发展计划项日作者简介:何以广(1984-),男,博土研究生,研究方向为高压状态方程;heyco6@mails,tsinghua,edu.cn1650强激光与粒子束第23卷13.1,e=10.8K,rn=0.29673nm,A=7796K,B=0.33467m-1,W=0.14mm。在经典分子动力学模拟过程中,每次模拟1000个原子,并采用微正则系综(NVE)。在每次模拟中,运行50000步,时间步长设为0.5fs。前25000步采用速度缩放稳定温度,从30000步开始热力学平均。2计算结果及讨论2.1氮的熔化曲线晶体熔化是一个复杂的转变过程。氦在高温高压区的熔化温度数据非常少,主要是因为静态高压技术不容易拓展到高压区。在低压区,采用经典分子动力学和经验势研究晶体的形态,利用 moldy软件模拟了1000个原子在等温等压(NPT)下的运动。在高压区,由于经验势不再适用采用了第一性原理分子动力A exp-6-MD学研究了氦的熔化温度。利用vasp软件,模拟了108…QMD个氦原子在2ps内的运动。对同一个压强,分别改变25温度,分析径向分布函数和分子扩散曲线,得到氦的熔化曲线,如图1所示。Fig. 1 Melting curve of heli在低压区,用exp-6势函数计算得到的熔化温度图1氮的熔化曲线比实验值低一些而用Aziz势函数计算得到的结果。与实验值1十分接近,这说明熔化温度对势函数非常敏感。在高压区,我们用第一性原理分子动力学得到的熔化温度与实验值在20GPa附近衔接得很好。在计算的范围内,氦的熔化温度随压强增大而增大。在50GPa时,氦的熔化温度Tm高达875K。值得注意的是,为了得到准确的熔化温度,采取的判据和采用的原子数目非常重要。在第一性原理分子动力学计算中,我们仅仅计算了108个氮原子,这可能导致一定的误差,为了得到更精确的结果需要采用更多的原子,但是这会导致计算量呈立方倍增加。 Delogu计算表明2,对铝来说,增大原子数,熔化温度有10%的变化。2.2氮的状态方程由于经验势函数来源于实验值,因此在低压区利用经典分子动力学和经验势来计算氦的状态方程是合适的。利用软件 moldy,我们计算了氮在温度T<0.86eV密度p<2.0g/cm2范围内的状态方程。在高压区,我们利用第一性原理分子动力学计算了氦的状态方程,利用软件vasp,计算得到了氦在温度T<4.5e、密度p<5.0g/cm2范围内的状态方程。图2是氦在低压区的 Hugoniot曲线、二次冲击波曲线2和T=0.3eV的等温曲线。 Hugoniot数据由Hugoniot方程得到E-E+0.5(p+p)(V-V)=0(4)可以看到,第一性原理分子动力学计算得到的 Hugonoit数据及二次冲击波数据在实验误差的范围内。图2中三角形和菱形数据点是氦T=0.3eV下的等温曲线,三角形是第一性原理分子动力学计算结果,菱形数据点FVT Ta2.73eV1000QMD second shockQMD, T=0.54ev40H-CMDT=0.3evMDT-03evdensity /(g-cm)Fig 2 Single- and double-shock Hugoniots of liquid heliumsFig 3 Isothe中国煤化工T=2.73visotherms of helium at T=0. 3 ev图3氦CNMHG等温曲线图2氯的 Hugoniot数据和二次冲击波数据及T=0.3eV等温曲线第6期何以广等:氦热力学性质的分子动力学研究1651是经典分子动力学计算结果。在低压区,二者十分接近,然而当压强在20GPa以上时,两种方法计算的结果出现偏差。图3比较了第一性原理分子动力学计算结果和化学模型FVT计算结果011.在低压区,二者接近,在高压区,二者出现较大的偏差。这可能是因为在高压区,FVT模型中各种近似不再成立。对计算结果进行了数值拟合得到了氮在温度T<4.5eV密度0.3g/cm2