DOI:10.19476/j.ysxb.1004.0609.2004.s1.058

LaNi₅-H体系计算机模拟

桑革 贾建平 沈崇雨 黄理 陈长安 武胜 孙颖

  中国工程物理研究院  

  中国工程物理研究院 绵阳621900  

摘 要：

采用CASTEP软件, 应用总能赝势方法对LaNi5与氢相互作用体系进行了模拟计算。计算中设置赝势为超软赝势, 对交换关联能项采用广义梯度近似 (GGA) , 对电子采用自旋极化处理。每个体系均采用完全的结构优化来计算平衡时的体系总能量。计算结果表明, 氢原子在钯晶格中最有利位置是八面体间隙位, 空位反而不利于氢原子占据。LaNi5H中氢原子在LaNi5晶格中最有利的位置是由2La2Ni组成的四面体间隙位。

关键词：

LaNi5;氢;模拟;

中图分类号： TP391.9

基金：中国工程物理研究院基金资助项目 (20020865);

Simulative calculations on LaNi₅-H

Abstract：

First-principles calculations on LaNi₅-H systems were made employing CASTEP software package. The ultrasoft pseudopotentials, generalized gradient approximation and electron spin polarized methods were used in calculations. A full geometry optimization was made to investigate the total energy of systems. From the calculated data it is found that the H atoms in vacancies are not stable although the volumes are bigger. H atoms are most prefer to diffuse in O-T-O pathway. The tetrahedral interstice consisting of 2La2Ni is the most stable site in perfect LaNi₅-H system. These two kinds of interfaces are all the most stable sites of H atoms. The H atoms will diffuse along the interface when trapped in it.

Keyword：

LaNi₅; hydrogen; simulation;

随着世界工业水平的日益提高, 人类对能源的需求也越来越大。 而传统的不可再生能源不但对自然环境存在严重的污染, 带来全球气候变暖等生态问题, 而且其储量有限, 已经很难满足人类下个世纪的能源需求。 所以, 寻找新的可再生且洁净的代替能源成了全人类迫在眉睫的任务。 研究氢能源利用技术的一个重要课题就是开发可以应用的储氢材料。 目前已有多种新材料纳入研究人员的视野。 最新的进展中有人提出以碳纳米管、 纳米纤维作为储氢材料, 也有很好的效果。 而传统的储氢金属因其特有的优点赢得了更为广泛的关注。 近年来, 有关其材料的制备或利用技术取得了长足的进展。

研究较多的另一种储氢材料是La-Ni合金体系。 La-Ni体系有多种化合物, 其中LaNi₅金属间化合物具有很强的吸氢能力。 LaNi₅是稀有储氢合金的典型代表, 它的优点是易活化、 平衡压力适中、 吸放氢动力学性能优异且滞后效应很小。 正因为有这些优点它被广泛应用于镍氢电池、 氢燃料电池, 能量转化与存储系统等。 LaNi₅的晶格结构为CaCu₅型, 在室温下一个晶格中能容纳6个氢原子, 生成具有六方晶格的LaNi₅H₆。 因而储氢量是比较可观的。

1 基本原理

CASTEP的基本原理仍然来自于量子力学。 众所周知, 材料在原子尺度上的行为遵循量子力学定律。 所以, 原则上可以用解Schrodinger方程的方法预测材料的性质。

1.1 密度泛函理论

密度泛函理论开始于Hohenberg-Kohn定理, 它提出系统所有的基态性质都是电荷密度的函数。 系统的总能量可以写成以下形式:

E_t[ρ]=T[ρ]+U[ρ]+E_xc[ρ] (1)

式中 T为无相互作用粒子的系统的动能; U为由于库仑相互作用的电子能; E_xc包括所有多体相互作用对总能的贡献, 这里特指交换和关联能。

本文通过波函数Ψ来建立电荷密度。 在分子轨道理论中, 波函数Ψ可以用单个粒子保函数来表示:

Ψ=A (n) |φ₁ (1) φ₂ (2) …φ_n (n) | (2)

这里的分子轨道是正交的:

〈φ_i|φ_j〉=δ_ij (3)

电荷密度可简单地由和的形式给出:

$\rho (r)={\displaystyle \underset{i}{\sum }|}\phi {}_i(r)|$ ² (4)

这里的波函数是所有被占据的分子轨道由方程 (4) 得出的密度即为电荷密度。

有了波函数和电荷密度, 则能量可写为

${\rm T}=?{\displaystyle \underset{i}{\overset{n}{\sum }}\phi }{}_i\left|\frac{-?{}^2}{2}\right|\phi {}_i$ 〉 (5)

$\begin{array}{l}U={\displaystyle \underset{i}{\overset{n}{\sum }}{\displaystyle \underset{\alpha }{\overset{{\rm N}}{\sum }}?}}\phi {}_i(r)\left|\frac{-{\rm Z}}{R{}_{\alpha }-r}\right|\phi {}_i(r)?+\\ \frac{1}{2}{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }?}}\phi {}_i(r{}_1)\phi {}_j(r{}_2)\frac{1}{r{}_1-r{}_2}\phi {}_i(r{}_1)\phi {}_j(r{}_2)?+\\ {\displaystyle \underset{\alpha }{\overset{{\rm N}}{\sum }}{\displaystyle \underset{\beta ＜\alpha }{\sum }\frac{{\rm Z}{}_{\alpha }{\rm Z}{}_{\beta }}{|R{}_{\alpha }-R{}_{\beta }|}}}.\\ \text{?}\text{?}=-{\displaystyle \underset{\alpha }{\overset{{\rm N}}{\sum }}?}\rho (r{}_1)\frac{{\rm Z}{}_{\alpha }}{|R{}_{\alpha }-r|}?+\\ \text{?}\frac{1}{2}?\rho (r{}_1)\rho (r{}_2)\frac{1}{|r{}_1-r{}_2|}?+\\ \text{?}{\displaystyle \underset{\alpha }{\overset{{\rm N}}{\sum }}{\displaystyle \underset{\alpha ＜\beta }{\sum }\frac{{\rm Z}{}_{\alpha }{\rm Z}{}_{\beta }}{|R{}_{\alpha }-R{}_{\beta }|}}}\\ =?-\rho (r{}_1)V{}_{{\rm N}}?+?\rho (r{}_1)\frac{V{}_{\text{e}}(r{}_1)}{2}?+V{}_{{\rm N}{\rm N}}\text{?}\text{?}\text{?}(6)\end{array}$

式中 Z_α指的是在N个原子的系统中核α所带的电荷。 等式右边第一项代表核与电子的相互吸引, 第二项代表电子间的相互排斥, 最后一项代表核与核之间的相互排斥。

式 (1) 中的最后一项即交换和关联能项在做计算时要做一定的近似。 一个简单却相当精确的近似是局域密度近似 (LDA) 。 局域密度近似假设在原子尺度上电子密度变化的相当缓慢 (每个分子区域就像是一个均一的电子气) 。 总的交换和关联能可以通过下面这个积分获得:

ε_xc[ρ]?∫ρ (r) ε_xc[ρ (r) ]dr (7)

式中 ε_xc[ρ]为均一电子气中每个粒子的交换关联能; r为粒子的编号。

考虑到电子气实际上并不均一, 则需要对局域自旋密度近似 (LSD) 模型进行改进。 这可以通过密度梯度展开来实现, 有时又被称作非局域自旋密度近似 (NLSD) 。 通常用到的有包括广义梯度近似 (GGA) 在内的很多种。

于是总能的表达式可以写为

$\begin{array}{l}E{}_{\text{t}}[\rho ]={\displaystyle \underset{i}{\sum }?}\phi {}_i\left|\frac{-?{}^2}{2}\right|\phi {}_i?+\\ \text{?}?\rho (r{}_1)[\epsilon {}_{\text{x}\text{c}}[\rho (r{}_1)]+\frac{V{}_{\text{e}}(r{}_1)}{2}-V{}_{{\rm N}}]?+V{}_{{\rm N}{\rm N}}\text{?}\text{?}\text{?}(8)\end{array}$

为了求解实际能量, 由式 (8) 可知:在正交化条件的限制下, 使式 (8) 取最小值的电荷密度ρ₀即为系统的基态电荷密度, 而此时的能量E_t即为系统的基态能量。

$\frac{\delta E{}_{\text{t}}}{\delta \rho }-{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }\epsilon }}{}_{ij}?\phi {}_i|\phi {}_j?$ =0 (9)

式中 ε_ij为拉格朗日乘数。

以上推导得出了一套Kohn-Sham方程:

$\left\{\frac{-?{}^2}{2}-V{}_{{\rm N}}+V{}_{\text{e}}+\mu {}_{\text{x}\text{c}}[\rho ]\right\}\phi {}_i=\epsilon {}_i\phi {}_i\text{?}\text{?}\text{?}(10)$

式中 μ_xc是交换关联势, 它是将E_xc微分后得到的。 对于局域自旋密度近似, 交换关联势可写作:

$\mu {}_{\text{x}\text{c}}=\frac{?}{?\rho }(\rho \epsilon {}_{\text{x}\text{c}})\text{?}\text{?}\text{?}(11)$

利用式 (10) 的本征值, 能量的表达式可变形为

$\begin{array}{l}E{}_{\text{t}}={\displaystyle \underset{i}{\sum }\epsilon }{}_i+?\rho (r{}_1)\{\begin{array}{l}\\ \epsilon {}_{\text{x}\text{c}}[\rho ]-\mu {}_{\text{x}\text{c}}[\rho ]-\\ \end{array}\\ \text{?}\text{?}\text{?}\frac{V{}_e(r{}_1)}{2}\}?+V{}_{{\rm N}{\rm N}}\text{?}\text{?}\text{?}(12)\end{array}$

实际应用中, 将分子轨道转化为原子轨道是比较方便的,

$\phi {}_i={\displaystyle \underset{u}{\sum }C}{}_{iu}\chi {}_u\text{?}\text{?}\text{?}(13)$

式中 原子轨道χ_u被称为原子基函数; C_iu为分子轨道展开系数。

基函数的选择多种多样, 有高斯函数、 Slater函数、 平面波、 数值轨道等。 其中CASTEP系统用的是平面波, 而Dmol³系统用的是数值轨道。

与分子轨道不同, 原子轨道并不正交。 这导致方程 (10) 变成了以下形式:

HC=εSC (14)

其中

${\rm H}{}_{uv}=?\chi {}_u(r{}_1)\left|\frac{-?{}^2}{2}-V{}_{{\rm N}}+V{}_{\text{e}}+\mu {}_{\text{x}\text{c}}\rho (r{}_1)\right|\chi {}_v(r{}_1)?\text{?}\text{?}\text{?}(15)$

S_uv=〈χ_u (r₁) |χ_v (r₁) 〉 (16)

1.2 超晶胞近似

超晶胞近似是CASTEP另一个重要的近似。

Bloch定理指出一个周期性系统中每一个电子的波函数都可以写成一个周期项和一个波函数项的乘积:

Ψ_i (r) =exp (ik·r) φ_i (r) (17)

式中 k为Bloch波矢; φ_i (r) 为包含晶体周期性的函数。

只有离散的一系列平面波拥有晶体的周期性, 所以, φ_i (r) 可以写成:

φ_i (r) =∑C_{i, G}exp (iG·R) (18)

于是, 通过Bloch定理每个电子状态都可以扩展为离散的平面波基函数组的形式。 它的主要好处是简化了Kohn-Sham方程:

${\displaystyle \sum [|k+G|{}^2\delta {}_{G{G}^{\prime }}+V{}_{\text{i}\text{o}\text{n}}(G-G^{\prime })+}$

$V{}_{\text{Η}}(G-G^{\prime })+V{}_{\text{x}\text{c}}(G-G^{\prime })]C{}_{i,k+G^{\prime }}=\epsilon {}_{i}G{}_{i,k+G}\text{?}\text{?}\text{?}(19)$

超晶胞近似的原则就是, 存在表面或缺陷的系统被看作是一个晶体的大的单元晶胞。

1.3 赝势

赝势的含义是考虑组合一系列在长度尺度上增加的量而忽略定义下的更长尺度下的相互作用。 在建立了赝势以后, 赝波函数即相当于核外价电子的波函数。 在离子芯区域内赝波函数是没有节点的, 相反真实的价电子波函数则必须有若干个节点。 所以赝波函数相对比较平滑。

LaNi₅是稀有储氢合金的典型代表, 它有易活化、 平衡压力适中、 吸放氢动力学性能优异且滞后效应很小的优点。 但是它也有自己的缺陷, 即易受到杂质气体的影响。 空气中的氧气或者少量的一氧化碳都能吸附在LaNi₅表面, 一部分O₂和CO直接与La反应生成氧化物和碳化物, 另一部分离解为O, C原子扩散进入合金。 表面氧化膜、 碳化膜阻止了合金进一步吸氢, 造成吸放氢速度变缓, 总的氢容量也会下降。 这种作用称为毒化作用。

对于这种新的储氢材料体系, 为了更为精确地设计制造实用的材料, 研究其微观性质是十分必要的。 然而目前的实验手段对微观领域的观测往往不尽人意, 所以作者考虑用CASTEP程序, 使用总能赝势方法从理论上计算整个结构的一些性质, 用以对材料的微观和宏观行为进行预测。

2 讨论

2.1 LaNi5单胞的能量和晶格常数

LaNi₅属于CaCu₅型六方晶系结构, 空间群为 $\text{D}{}_{6\text{h}}^1-\text{Ρ}\frac{6}{\text{m}}$ mm, 其结构如图1所示。

图1 LaNi5单胞的空间构型

每个单胞中有6个原子, 1个La原子占据1 (a) 位置, 原子坐标为 (0, 0, 0) , 2个Ni原子 (A位) 占据2 (c) 位置, 原子坐标为 (1/3, 2/3, 0) 和 (2/3, 1/3, 0) , 3个Ni原子 (B位) 占据3 (g) 位置, 原子坐标为 (1/2, 0, 1/2) 、 (0, 1/2, 1/2) 和 (1/2, 1/2, 1/2) 。 点阵常数的实验值为a=0.501 7 nm, c=0.398 1 nm, V= 8.678×10^-2 nm³。 用该实验值带入程序构建LaNi₅晶胞, 然后降低其对称性为P1。 采用CASTEP程序进行结构优化。 考虑到LaNi₅晶胞所含原子数较多, 结构较为复杂, 为减少计算所需时间, 设置运算精度为coarse, 即能量变化小于5×10^-5时即可认为是收敛的。 优化的结果如下:

E_t=-7 648.184 2 eV,

a=0.508 12 nm, b=0.505 33 nm,

c=0.435 22 nm;

α=90.055 3°, β=90.029 7°, γ=120.772 2°;

V=9.601 61×10^-2 nm³

2.2 氢原子在间隙位置的能量

LaNi₅晶格中可认为存在3种间隙位, 图2所示分别为2La2Ni构成的四面体间隙位 (T₁) , 4Ni构成的四面体间隙位 (T₂) 和2La4Ni构成的八面体间隙位 (O) 。

图2 LaNi5单胞中的间隙位

在这3种间隙位中分别引入1个氢原子, 即相当于LaNi₅H的结构, 进行构型优化。 其晶胞分别见图3～5。优化结果列于表1。

由表1中含间隙氢原子的单胞总能量数据减去完好LaNi₅单胞的总能量即为氢原子在晶格间隙位置的能量, 它与自由氢原子的能量差即为形成能ΔE ${}_{\text{Ι}}^{\text{f}}$ , 即

ΔE ${}_{\text{Ι}}^{\text{f}}$ =E ${}_{\text{Η},\text{Ι}}^{\text{t}}$ -E ${}_{\text{Η},\text{f}\text{r}\text{e}\text{e}}^{\text{t}}$ =E ${}_{\text{L}\text{a}\text{Ν}\text{i}{}_5,\text{Η}}^{\text{t}}$ -E ${}_{\text{L}\text{a}\text{Ν}\text{i}{}_5}^{\text{t}}$ -E ${}_{\text{Η},\text{f}\text{r}\text{e}\text{e}}^{\text{t}}$ (24)

式中 E ${}_{\text{Η},\text{Ι}}^{\text{t}}$ , E ${}_{\text{Η},\text{f}\text{r}\text{e}\text{e}}^{\text{t}}$ , E ${}_{\text{L}\text{a}\text{Ν}\text{i}{}_5,\text{Η}}^{\text{t}}$ , E ${}_{\text{L}\text{a}\text{Ν}\text{i}{}_5}^{\text{t}}$ 分别为间隙氢原子、 自由氢原子、 含氢原子的LaNi₅单胞和完好LaNi₅单胞的总能量值。 所以由表1中的数据和上一章计算得出的自由氢原子能量数据, 可以得出各间隙氢原子的形成能分别为

图3 T1位含1个氢原子的LaNi5晶胞图像

图4 T2位含1个氢原子的LaNi5晶胞图像

图5 O位含1个氢原子的LaNi5晶胞图像

ΔE ${}_{\text{Ι}}^{\text{f}}$ (四面体T₁) =-2.965 4 eV,

ΔE ${}_{\text{Ι}}^{\text{f}}$ (四面体T₂) =-2.461 4 eV,

ΔE ${}_{\text{Ι}}^{\text{f}}$ (八面体O) =-2.397 8 eV。

1个LaNi₅晶格中添加1个氢原子的理论生成热可由下式计算:

$\begin{array}{l}\text{Δ}{\rm H}=E{}_{\text{Η},\text{Ι}}^{\text{t}}-\frac{1}{2}E{}_{\text{Η}{}_2,\text{f}\text{r}\text{e}\text{e}}^{\text{t}}=\\ E{}_{\text{L}\text{a}\text{Ν}\text{i}{}_5,\text{Η}}^{\text{t}}-E{}_{\text{L}\text{a}\text{Ν}\text{i}{}_5}^{\text{t}}-\frac{1}{2}E{}_{\text{Η}{}_2,\text{f}\text{r}\text{e}\text{e}}^{\text{t}}\text{?}\text{?}\text{?}(25)\end{array}$

表1 LaNi5晶胞中添加1个间隙氢原子后的总能量和晶格常数的优化结果

晶胞	总能量/eV	a/nm	b/nm	c/nm	α/ (°)	β/ (°)	γ/ (°)	V/nm3
初始单胞	-7 648.184 2	0.508 12	0.505 33	0.435 22	90.055 3	90.029 7	120.772 2	9.601 61×10-2
间隙T1	-7 664.720 7	0.544 31	0.539 83	0.428 13	89.744 4	90.280 1	123.901 2	10.441 36×10-2
间隙T2	-7 664.216 7	0.535 02	0.532 08	0.430 74	90.439 0	89.787 6	122.191 5	10.376 66×10-2
间隙O	-7 664.153 1	0.529 45	0.534 97	0.422 20	89.874 1	90.024 5	121.517 5	10.194 33×10-2

式中 E ${}_{\text{Η}{}_2,\text{f}\text{r}\text{e}\text{e}}^{\text{t}}$ 为自由氢分子的总能量。

由表1中的数据和本文计算得出的自由氢分子能量数据, 可以得出各间隙氢原子的理论生成热分别为

ΔH (四面体T₁) =-0.714 4 eV,

ΔH (四面体T₁) =-0.210 4 eV,

ΔH (八面体O) =-0.146 8 eV。

另外, 由表1中数据可见, 在单胞中插入杂质原子后, 晶格点阵会发生一定的畸变, 体积膨胀。 模拟结果显示单胞体积膨胀分别为

ΔV/V (四面体T₁) =8.7%,

ΔV/V (四面体T₁) =8.1%,

ΔV/V (八面体O) =6.2%。

纯以能量的观点来考虑, 则根据模拟结果, 氢原子在LaNi₅单胞中最可能的位置是四面体T₁位, 其次是四面体T₂位, 最后是八面体位。

3 结论

1) 基于密度泛函理论、 广义梯度近似、 赝势平面波方法, 可以计算出单个氢原子在LaNi₅晶格各间隙位中的形成能数据。

2) 在整个晶体范围内最有利于容纳氢原子的是由2个La原子和2个Ni原子构成的四面体间隙位T₁, 其次是由4个Ni原子构成的四面体间隙位T₂, 最后是由2个La原子和4个Ni原子构成的八面体间隙位O。

参考文献

[1]　DavisWD.KnollsAtomicPowerLabRep, 1954, 1227.

[2]　WorshamJE, WilkisonJRMK, ShullCG.JPhysChemSolidsPergamonPress, 1957 (3) :303310.

[3]　BohmholdtG, WicheE.ZPhysChemNF, 1967, 56:133.

[4]　GillanMJ.JPhysC, 1990, 19:6169.

[5]　SalomonsE.JPhysCondensMatter, 1990, 2:845.

[6]　陈立新.中国有色金属学报, 1999, 9 (1) :61.

[7]　MilmanV, PayneMC.PhysicalReviewLetter, 1993, 10.

[8]　YamagishiS, JenkinsSJ.JournalofChemicalPhysics, 114 (13) :

[9]　PhilipJD.LindanNM.Harrison, PhysicalReviewLetters, 1998, 26:

[10]　PughS, GillanMJ.SurfaceScience, 1994, 320:331343.

[11]　GoniakowskiJ, GillanMJ.SurfaceScience, 1996, 350:145158.

[12]　SzymanskiMA, GillanMJ.SurfaceScience, 1996, 367:135148.

[13]　KantororichLN, GillanMJ.SurfaceScience, 1997, 374:373386.

[14]　KuaJ.JournalofChemicalPhysics, 2001, 22