0 引言

近年来，针对具有A₅B₃型化合物的研究越来越多，这是因为该类化合物具有独特的物理化学性能^[1-5]，研究表明超过40个金属间化合物具有A₅B₃型结构。常见的结构类型有Cr₅B₃型、W₅Si₃型和Mn₅Si₃型这3类，Tao等^[5-8]利用基于密度泛函的第一性原理方法对Ta₅Si₃、W₅Si₃、V₅Si₃以及Tl₅Te₃等化合物的物理化学性能进行模拟计算研究，获得了这些化合物的结构稳定性、力学性能、电子结构特性以及热物理性能和热电性能等。含碱土金属Ca的化合物中，Ca₅Si₃因其有潜在的储氢性能^[3]和热电性能^[9]备受关注。然而在含Ca的化合物中，Ca₅Zn₃也具有Cr₅B₃型结构^[1]。2012年，Yang等^[10]利用基于密度泛函的第一性原理的方法对Ca-Zn体系的7个化合物进行了计算研究，计算获得了形成焓、弹性常数和电子态密度等信息。然而，至今还没有针对Ca₅Zn₃的更加系统研究的报道。为了能够获得更加系统的物理性能，本文将利用基于密度泛函的第一性原理方法对Ca₅Zn₃化合物进行系统计算研究，获得该化合物的热力学稳定性、力学性能、电子结构特性，并结合德拜模型对其热物理性能进行研究和预测。

1 材料与方法 1.1 Ca₅Zn₃的晶体结构

采用Cr₅B₃型、W₅Si₃型和Mn₅Si₃型等3种结构来讨论Ca₅Zn₃结构稳定性。其中Cr₅B₃型结构是一个四方晶系的结构，Pearson符号为tI32，共有4种不同占位，Cr占据16l和4c位，B占据4a和8h位。每个Cr₅B₃型结构单胞中含有10个Cr原子，6个B原子。其余两种结构的具体信息可参见文献[5]。3种晶体结构如图 1所示。

1.2 计算方法

本文采用基于密度泛函理论(Density Functional Theory，DFT)的第一性原理计算软件包Vienna Ab initio Simulation Package (VASP)^[11]进行计算。第一性原理计算方法是投影缀加波方法^[12-13]，交换关联势采用广义梯度近似的Perdew-Burke-Ernzerhof (PBE)方案^[14]。布里渊区积分采用Monkhorst Pack布点方法^[15]，对于Cr₅B₃和W₅Si₃结构，K点取11×11×11，而对Mn₅Si₃结构，K点为9×9×11。在本文计算中, 波函数的截断能设置为500 eV，能量收敛精度为10^-6 eV。在进行性能计算前，所有结构都进行全优化，直到原子间的力小于0.01 eV/Å为止，以保证计算的准确性。有关形成焓、弹性常数以及热物理性质的计算方法参见文献[5]。

Ca₅Zn₃化合物形成焓可以用以下的公式进行计算：

$ \begin{array}{c} \mathit{\Delta} H\left(\mathrm{Ca}_{5} \mathrm{Zn}_{3}\right)=E_{\text {total }}\left(\mathrm{Ca}_{5} \mathrm{Zn}_{3}\right)-5 E_{\text {total }}(\mathrm{Ca}) \\ 3 E_{\text {total }}(\mathrm{Zn}), \end{array} $

(1)

其中，E_total(Ca₅Zn₃)、E_total(Ca)和E_total(Zn)分别是Ca₅Zn₃合金、Ca以及Zn元素的总能量。

晶体弹性常数是一种反映材料基本性质的常数，它描述了晶体对外加应变的响应刚度。当忽略温度对体系总能的影响时，在应变很小的情况下，体系的内能与应变的大小存在二次线性关系(胡克定律)，弹性常数就是描述这种二次线性关系，即二次线性项的系数。为了计算弹性常数，我们通过设定一组形变量，求解固体在不同形变下所对应的总能量，利用弹性理论公式拟合这组能量就可以得到它们的弹性常数。通过定义一个应变张量ε，可使基矢$ \overrightarrow a $转换为新矢量：

$ {{\vec a}^\prime } = (\vec I + \varepsilon )\vec a, $

(2)

$ \varepsilon = \left( {\begin{array}{*{20}{c}} {{\varepsilon _1}}&{{\varepsilon _6}/2}&{{\varepsilon _5}/2}\\ {{\varepsilon _6}/2}&{{\varepsilon _2}}&{{\varepsilon _4}/2}\\ {{\varepsilon _5}/2}&{{\varepsilon _4}/2}&{{\varepsilon _2}} \end{array}} \right), $

(3)

$ \overrightarrow I $为3×3的单位矢量。采用Voigt标记xx→1、yy→2、zz→3、yz→4、zx→5、xy→6。在这套符号系统中，弹性常数被写为C_ij，i和j遍及1—6。定义应力张量σ_i，由此，体系总能的Taylor展开可以表示为

$ \begin{array}{*{20}{c}} {E(V,\varepsilon ) = E\left( {{V_0},0} \right) + {V_0}\sum\limits_i {{\sigma _i}} {\varepsilon _i} + }\\ {\frac{{{V_0}}}{2}\sum\limits_{ij} {{C_{ij}}} {\varepsilon _i}{\varepsilon _j} + O\left( {\varepsilon _i^3} \right),} \end{array} $

(4)

其中，V₀是晶格平衡时的体积。上式第二项是一个线性项，如果体积没有变化或者在平衡位置处，该项为零；O(ε_i³)是在展开的多项式中忽略ε_i的立方项或者更高次项。

Cr₅B₃结构类型的Ca₅Zn₃是四方晶系。四方晶系有6个独立弹性常数，分别为C₁₁, C₁₂, C₁₃, C₃₃, C₄₄和C₆₆。通过下面6个形变矩阵可以求得相应弹性常数(表 1)。

德拜温度作为一个基本的物理参数，与热容、弹性常数等有着密切联系。一般来说，德拜温度越高，化合物熔点越高，意味着原子间的相互作用越强。一般来说，德拜温度(Θ)可以从平均声速(u_m)进行估算:

$ \mathit{\Theta}=\frac{h}{k_{B}}\left(\frac{3}{4 \pi V}\right)^{1 / 3} u_{m}, $

(5)

其中, h和k_B分别是普朗克常数和玻尔兹曼常数，V是原子平均体积，在多晶体中的平均声速(u_m)可以用以下公式计算:

$ u_{m}=\left[\frac{1}{3}\left(\frac{2}{u_{t}^{3}}+\frac{1}{u_{l}^{3}}\right)\right]^{-1 / 3}，$

(6)

其中，u_t和u_l分别是横向声速和纵向声速，这两种声速可以通过以下公式获得:

$ {{u_t} = {{\left( {\frac{G}{\rho }} \right)}^{1/2}},} $

(7)

$ {{u_l} = {{\left( {\frac{{3B + 4G}}{{3\rho }}} \right)}^{1/2}},} $

(8)

其中, G为剪切模量，B为体积模量，ρ为质量密度。

2 结果与分析 2.1 结构稳定性能

为了获得Ca₅Zn₃的结构稳定性，考虑了Cr₅B₃型、W₅Si₃型和Mn₅Si₃型3种晶体结构，计算获得的晶格常数、形成焓(表 2)，同时分析3种晶体结构的总能随体积变化的关系(图 2)。本文计算获得的Cr₅B₃型晶格常数为7.865 6 Å和15.366 3 Å，Bottcher等^[1]的实验值为7.963 Å和15.407 Å，其他第一性原理计算值^[10]为7.802 Å和15.620 Å。本文的晶格常数计算值与实验值之间的误差为1.22%和0.26%，比Yang等^[10]计算值的误差(2.02%，1.38%)小，说明本文计算获得的晶格常数更为合理。W₅Si₃型和Mn₅Si₃型Ca₅Zn₃化合物的晶格常数也作为预测结果列于表 2。从表 2中给出的形成焓可以看出，Cr₅B₃型Ca₅Zn₃的形成焓最低为－0.201 eV/atom，其他两种晶体结构的形成焓分为－0.160 eV/atom和－0.164 eV/atom。从能量的角度，其他两种晶体结构的形成焓比Cr₅B₃型的形成焓高出0.04 eV/atom，因此，具有Cr₅B₃结构的Ca₅Zn₃是最稳定的化合物。从图 2中可以看出，具有Cr₅B₃型的Ca₅Zn₃具有最低的总能量，同样表明Cr₅B₃型Ca₅Zn₃是最稳定的化合物，这与Bottcher等^[1]的实验结果是一致的。

2.2 力学性能

为了研究Ca₅Zn₃的力学性能，本文计算了稳定相的晶格常数，由上文可知Cr₅B₃型晶体结构属于四方晶系，其独立的弹性常数为6个，为此，本文利用应力应变的方法^[16]对其弹性常数进行了计算(表 3)。根据四方晶系力学稳定性条件：C₁₁>0, C₃₃>0, C₄₄>0, C₆₆>0, (C₁₁－C₁₂)>0, (C₁₁+C₃₃－2C₁₃)>0和(C₁₁+C₁₂)C₃₃－2C₁₃²>0，本文计算获得的Ca₅Zn₃的弹性常数都满足力学稳定性条件，表明该化合物具有静力学稳定性。C₁₁大于C₃₃，由此可见在[100]和[010]方向上的键强度高于[001]方向上的键强度。多晶体的体积模量、剪切模量以及杨氏模量都可以从弹性常数给出^[5]。本文计算得到的体积模量和Yang等^[10]报道的计算值吻合，而剪切模量的计算值比Yang等^[10]报道的计算值稍小，这是由于本文计算的C₃₃比Yang等^[10]报道的计算值稍小导致的。

B/G值一般用来预测化合物的脆韧性^[17]，该比值大于1.75，化合物呈现韧性，比值越大韧性越大；反之，该比值小于1.75，则化合物呈现脆性，且比值越小脆性越大。Ca₅Zn₃的B/G值为2.61，说明该化合物呈现韧性。为了研究化合物的各向异性，对于四方晶系来说，一般用A₁=2C₆₆/(C₁₁－C₁₂)和A₂=2C₄₄/(C₁₁－C₁₂)进行说明。该数值为1时，则呈现各向同性，偏离1越远则各向异性越明显。为了更加明显地呈现该化合物的各向异性，图 3给出了杨氏模量的空间分布图，由此可以清晰地看出在不同的方向上，杨氏模量的数值存在明显的差异。

2.3 电子结构特性

为了深入了解该化合物原子间的相互作用，本文计算了该化合物的电子态密度和电荷密度差(图 4)。由图 4a可以看出，在费米面以下的电子态密度，主要是Ca的s电子和p电子以及Zn的s电子和p电子的贡献，费米面以上的电子态密度主要是Ca的d电子和Zn的p电子的贡献。在费米面附近，Ca的p电子和d电子跟Zn的p电子发生杂化效应。在费米能外态密度不为零，说明该化合物呈现金属性。由图 4b和图 4c可知，Ca和Zn之间主要呈现离子键特性，Ca失去电子，而Zn则获得电子，这个和图 4a中的电子态密度在费米面附近发生杂化效应是一致的。由于Ca的电负性为1.0，Zn的电负性为1.6，电负性越大表明越容易获得电子，由此可见本文的计算结果是合理的。

2.4 热物理性能

Ca₅Zn₃化合物的德拜温度为202 K。平均声速为2 097 m/s，纵向声速为3 072 m/s，横向声速为1 901 m/s (表 3)。利用德拜模型^[18-21]，结合第一性原理总能量的计算结果，计算了Ca₅Zn₃化合物在0—1 000 K温度范围、0—10 GPa压强范围内的热物理性质(图 5)。由图 5a可知，化合物的体积随着温度的升高而增加，但是随着压强的增加则降低。并且随着压强的越来越大，体积减小的幅度越来越小，说明原子间的作用力也越来越大。图 5b给出了体积模量随温度和压强的变化关系，随着温度的升高，体积模量有所下降，但是随着压强的增加，体积模量也在增加。图 5c表示热膨胀系数随着温度的升高而增加，随着压强的增大而热膨胀系数增加的幅度降低。在6 GPa以上时，随着温度的上升，热膨胀系数增加的幅度极小。表明在高压下，热膨胀系数变化小，体积膨胀不明显。图 5d展示的是定容热容随温度和压强的变化关系，在0—300 K区间，定容热容急速增加，在300—400 K区间，等容热容的上升趋势变缓，在400 K以上，已经接近等容热容的极限值，即杜隆珀蒂极限。

3 结论

1) 具有Cr₅B₃型的Ca₅Zn₃形成焓最低，是最稳定的结构，与前人的实验结果一致。

2) 弹性常数的计算值显示其满足四方晶系的力学稳定性标准，计算得到了体积模量、剪切模量、杨氏模量以及B/G值，B/G值显示该化合物呈现韧性特征。杨氏模量呈现各向异性。

3) 电子态密度显示Ca的p电子和d电子与Zn的p电子在费米能附近发生杂化效应，费米面上有态密度说明该化合物为金属性。电荷密度差显示Ca失去电子、Zn得到电子，与这两种元素的电负性符合，同时呈现离子键特性。

4) 体积、体积模量、热膨胀系数以及定容热容随着温度和压强的变化趋势说明，随着压强的增加、体积变小，体积模量变大、热膨胀系数变小，定容热容在低温下变小，高温下趋于一致。