材料的非简谐性描述符*

李环娅 周柯 尹万健

(苏州大学能源学院,能源与材料创新研究院,江苏省先进碳材料与可穿戴能源技术重点实验室,苏州 215006)

非简谐效应是诸如软模相变、负热膨胀、多铁性和超低热导率等材料性质的根源.已有的关于量化材料非简谐性的方法没有给出清晰准确的非简谐性描述符,并且计算流程复杂,需要极其耗时的分子动力学模拟.故亟需提出一个可以快速计算的非简谐性描述符,用来理解、评估、设计和筛选具有强非简谐性的功能材料.本研究将晶格非谐性分解为单声子非谐性 ,并提出温度依赖的晶格非谐性的定量描述符 Aph(T) .该描述符既可以定量描述从Si,GaAs,CdTe,NaCl 到CsPbI3 的晶格非谐性的变化趋势,又可以成功预测非简谐效应驱动的体积模量和晶格热导率性质的变化.本工作提出的非简谐性描述符能够快速量化材料非简谐性,并且可直观地展现材料非简谐效应的声子模态分布.本方法计算简单、高效且有效,可为基于非简谐性筛选与设计材料打下基础.

1 引言

近年来,凝聚态体系涌现了诸多具有特定性质的功能材料,如负热膨胀材料[1-4]、铁电材料[5,6]、超低热导率材料[7-11]、光电材料[12-14]等.这些材料具有的特殊性质往往起源于非简谐效应,无法通过简谐理论来解释.以其中的光电功能材料为例,对于传统的光电材料如Si,GaAs,CdTe 等,因其非简谐性较弱,所以以往的研究大部分着眼于其晶体缺陷,能带等角度来解释影响其光电转换效率的主要因素[15-17],而并没有关注其非简谐效应.但新一代光电材料卤化物钙钛矿就是一种典型的强非简谐功能材料,不同于传统的半导体,如Si,GaAs和CdTe 等具有刚性的共价网络,卤化物钙钛矿的晶格非常“软”[18-20].最近一些研究表明,钙钛矿的软晶格可能是导致其独特的载流子传输机制、长的载流子寿命和独特的离子迁移率机制的原因[21,22].钙钛矿晶体中存在低频、不稳定、八面体倾斜模式的无序振动,其A 位阳离子的振动模式和无机BX3-骨架的振动耦合导致其具有极强的非简谐性,进而导致其软晶格的形成[7,23-28].有研究认为,钙钛矿软晶格的非简谐振动抑制了其载流子的离域,使载流子自陷成大的极化子,导致其复合的能垒增大,进而降低载流子复合速率,这是钙钛矿材料具有长载流子寿命特性的主要原因[29-31].钙钛矿的非简谐性也会导致钙钛矿材料中B—X 键的减弱,从而加速离子迁移[32-36].对于像钙钛矿这种强非简谐性的光电功能材料,对其研究不能如以往的对Si,GaAs,CdTe 等弱非简谐材料一样忽略非简谐性对其功能性质的影响.准确评估材料的非简谐性,不仅有助于研究人员理解、评估其非简谐性的作用机理,甚至可以帮助其设计和筛选优异的强非简谐太阳能电池材料.

目前材料的非简谐性常用第一性原理结合分子动力学(molecular dynamics,MD)来进行计算模拟,包括自洽声子理论(Self-consistent phonon theory,SCP)[37-40]、自洽场理论(Self-consistent field,SCF)[41,42]、基于从头算分子动力学(ab-initiomolecular dynamics,AIMD)的变温有效势(temperature dependent effective potential,TDEP)[7,43,44]等.模拟结果常由非抛物线冻结声子势、体系自由能修正、三阶原子间力常数、声子频移和弛豫图来表示.然而,这些方法没有明确地提出一个非简谐性描述符,并且需要繁杂和耗时的MD 计算.最近,Knoop等[13]提出了一种表征非简谐性的新方法,该方法虽然创新地提出了一个以描述对象受力为基础的非简谐性描述符,但其仍然依赖于复杂耗时的AIMD 计算.

因为材料的特定性质往往只与该体系中特定的声子振动模式相关,为了更好地理解非简谐效应对该特定性质的影响,探寻其作用机理,本工作选取声子作为衡量非简谐性强弱的研究对象,以便单独分析影响某个特定性质的特定振动模式.本工作中推导并讨论了一种基于声子理论和静态计算的统计方法,该方法通过计算一个新的非简谐性描述符来量化材料的晶格非简谐性.该方法将晶格非简谐性分解为独立声子振动模式的非简谐性贡献,提出了声子的定量非简谐性描述符,并将其与玻色-爱因斯坦分布结合起来提出了温度依赖的材料非简谐性的定量描述符Aph(T) .选取具有代表性的单质和化合物体系作为研究对象,以材料体模量和室温下晶格热导率为参考,验证了该描述符的正确性.该描述符的计算过程不需要进行昂贵的分子动力学计算,能够快速量化材料的非简谐性,从而对如光电材料一样受非简谐性影响而具有特定性质的功能材料的高通量分类、设计和筛选提供可能.

2 计算模型和方法

2.1 声子非简谐性描述符

体系的声子振动对自由能的贡献E可以分解为单个声子振动贡献的总和:

式中V(Rq,j) 是单个声子的振动对自由能的贡献,Rq,j表示该声子振动模式中原胞原子偏离其平衡位置的集体位移.为了消除量纲影响,首先对振动模式Rq,j进行归一化处理,即αq,j=Rq,j/|Rq,j|,令αq,j∈[0,1](为无量纲数),然后手动设定rq,j=d·αq,j,d为我们所设定的原子位移的偏离量,当所有原子都按原振动模式中的位移比例乘以d偏移时,V(rq,j) 可以用泰勒展开式表示为

当d取统一值时,∅n只与声子的振动模式有关,其值不随温度的变化而变化,是声子的本征属性.因此,可定义为声子的非简谐描述符:

图1 为钙钛矿原胞中300 K 下两种能够被激发的声子振动模式的非简谐性描述符的示意图 (虚频声子的示意图详见补充材料(online)中图S1).图1(a)中红点和蓝点的符合度较高,其声子的非简谐性描述符小于图1(b),即图1(a)中声子振动模式的非简谐性弱于图1(b)中的声子振动模式.

图1 不同振动模式的声子的非简谐性示意图,图中红色圆点表示完全简谐情况下声子振动对体系自由能的贡献 V (rq,j) 随着振动幅度 |rq,j| 的变化趋势,蓝色圆点则表示真实情况下其变化趋势Fig.1.The schematic diagram of phonon anharmonicity of different vibration modes.The red dots represent the change trend of the contribution of phonon vibration to the free energy of the system V (rq,j) with the vibration amplitude |rq,j| in harmonic approximation,and the blue dots represent the change trend in the real case.

2.2 材料非简谐性描述符

声子为玻色子,其在热平衡状态下服从玻色-爱因斯坦统计.体系中不同振动模式的声子占据数与其频率关系为

其中ωq,j代表声子的频率,T代表温度.因此,理论上体系中所有被激发声子的非简谐贡献的加和即为该温度下该体系的总非简谐贡献.在此定义Aph(T)为材料的非简谐性描述符:

Aph(T) 是依赖于温度的材料的本征属性,N是体系中声子振动模式的总取样数.

以Ge 元素为中心,挑选了弱非简谐性的C,Si,Ge,Sn,GaAs,CdTe,NaCl 晶体(如图2(a)所示)和强非简谐性的β-CsPbI3和γ-CsPbI3晶体作为研究对象,其组成元素在元素周期表中横纵向均涵盖了较广的范围.所选晶体的结构如图2(b)所示,C,Si,Ge,Sn 的空间群均为,而GaAs,CdTe,NaCl 的空间群均为.这7 种晶体用于验证我们所提出的非简谐性描述符对于量化同构型不同元素组成的晶体的非简谐性的适用性.β-CsPbI3的空间群为I4/mcm,γ-CsPbI3的空间群为Pnma,这两者用于验证我们所提出的非简谐性描述符对于量化同元素组成不同构型的晶体的非简谐性的适用性.

图2 (a) 所选材料的组成元素分布;(b) Fd¯3m (C,Si,Ge,Sn);(c) F¯43m (GaAs,CdTe,NaCl);(d) I4/mcm ((β-CsPbI3);(e) Pnma(γ-CsPbI3);(f)计算流程图Fig.2.(a) The distribution of constituent elements of typical materials we selected;(b) Fd¯3m (C,Si,Ge,Sn);(c) F¯43m (GaAs,CdTe,NaCl);(d) I4/mcm (β-CsPbI3);(e) Pnma (γ-CsPbI3);(f) calculation flow diagram.

计算流程如图2(c)所示,首先利用Phonopy软件计算出所选晶体的声子谱,然后取样分析单个声子的振动模式Rq,j对应的原胞中原子的位移矢量,并将归一化后的位移矢量比例乘以位移d得到新的位移矢量,然后将晶体原胞中的原子按照新的位移矢量手动调整位置,直至原胞中新位移矢量中的最大位移偏离达到我们所设定的偏离标准.因300 K 下原子振动位移一般不超过键长的1/10,所以将偏离标准设置为0.3 Å,这个大小在经验上既不会过大到对真实情况的模拟失真,也不会过小到忽略非简谐效应.随后我们通过对手动调整偏离后的原胞进行静态自洽计算来得到该体系的能量,就可以得到不同d值对应的体系能量,即该声子模式的振动对自由能的贡献.可以通过最小二乘法得到该声子的振动能量-位移曲线.最后,由(4)式得到单个声子的非简谐描述符,由(6)式得到与温度相关的材料非简谐描述符Aph(T) .

本工作中密度泛函理论(density functional theory,DFT)计算均采用VASP 软件包实现.我们使用投影缀加波 (PAW)和Perdew-Burke-Ernzerhof (PBE)交换关联泛函来对核与电子之间的相互作用进行模拟.其中结构优化和静态计算所用的平面波截止能量为500 eV,能量收敛精度则分别设为10-6eV 和10-8eV,C,Si,Ge,Sn,GaAs,CdTe 和NaCl 晶体对应的布里渊区k点网格选取为5×5×5,γ-CsPbI3和β-CsPbl3晶体对应的布里渊区k点网格选取分别为3×2×3 和3×3×4.其中声子谱采用密度泛函微扰理论(density functional perturbation theory,DFPT)方法计算,其计算采用Phonopy 软件[45,46].晶格热导率计算则采用Phono3py[46,47].其中声子谱和热导率计算所用的平面波截断能为500 eV,能量收敛精度为10-8eV,C,Si,Ge,Sn,GaAs,CdTe 和NaCl 晶体的超胞尺寸分别为2×2×2,原子总数为16,对应的布里渊区k点网格选取为5×5×5,γ-CsPbI3和β-CsPbl3晶体的超胞尺寸为1×1×1,原子总数分别为20 和10,对应的布里渊区k点网格选取分别为3×2×3 和3×3×4.其中体模量计算中所用的平面波截止能量为500 eV,能量收敛精度则为10-5eV,施加的弹性应变小于2%,选取的一组应变值分别为-2%,-1%,0%,1%,2%.

3 结果和讨论

我们计算了所选9 种材料的声子谱,然后在完整的布里渊区路径上等间隔取样计算该点声子模态的非简谐性描述符,并将计算出的样点声子的值作为颜色值标注在声子谱上.为方便展示,将横轴截取为Γ 点到原布里渊区路径的前1/4,纵轴截取到4.3 THz (300 K 下能被激发的声子的最高频率)然后汇总成图3 (完整的声子谱非简谐模态分布图详见补充材料(online)图S2—图S10).通过对数据的分析,研究发现图3 中存在的简并声子模态的非简谐性描述符值相同.如图3 所示,该描述符可以直观地展示材料非简谐性在整个声子谱中不同模态的分布,从而为实现与特定性质相关的声子模态的独立分析提供可能.从图3 可以看出,当声子的频率分布接近时,例如图3(e),(g),可以直接从图中点的颜色判断该材料非简谐性的强弱.而当点的颜色相近时,例如图3(e),(f),其声子的比较接近,这时就需要考虑声子的频率分布对总非简谐性的影响.

图3 所选材料的声子谱,图中点的颜色值为对应声子模态的非简谐性描述符 (a) C;(b) Si;(c) Ge;(d) Sn;(e) GaAs;(f) CdTe;(g) NaCl;(h) γ-CsPbI3;(i) β-CsPbI3Fig.3.The phonon spectrum of the selected material,the color value of the point in the graph is the anharmonic descriptor of the corresponding phonon mode:(a) C;(b) Si;(c) Ge;(d) Sn;(e) GaAs;(f) CdTe;(g) NaCl;(h) γ-CsPbI3;(i) β-CsPbI3.

图4(a)展示了GaAs 和CdTe 的一条横模声学声子支(TA)的的差值,从图3(e),(f)可以看出,该声子支在两种材料中的颜色接近,从图4(a)也可以得出两者数值差异不大,但是从图3(f)可以看出,CdTe 中低频声子占据数远高于图3(e)中的GaAs,因此总非简谐性CdTe 强于GaAs.图4(b)统计了我们所选材料在Γ 点附近的低频声子的非简谐性贡献和所有路径上的声子非简谐贡献总和.如图4(b)所示,对9 种材料中Γ 点附近低频声子的非简谐贡献的分析显示,随着非简谐性的增加,Γ 点附近的低频声子对整体非谐性贡献占比逐渐减小.因此结合图4(a),(b),可以得出材料的非简谐性既取决于声子的非简谐性,也取决于给定温度下被激发的声子占比分布.我们提出的描述符Aph(T) 同时考虑了这两个因素,为评估温度依赖的材料的非简谐性提供了一个清晰而准确的标准.

目前一些研究证明,非简谐性往往是材料软晶格形成的原因[51-53],晶格的软硬程度则可以用体模量来进行衡量.具有较强非简谐性的材料往往具有较小的体模量,弱非简谐材料则具有较大的体模量[53-56].因此我们分别计算了室温下所选9 种材料的晶体非简谐性描述符和体模量来验证此描述符的正确性.C,Si,Ge,Sn 是以Ge 为中心,并沿周期表纵向延伸的同族元素组成的材料,随着晶体的非简谐性逐渐增大,这4 种材料的体模量逐渐减小.我们对其的计算表明,Aph(300 K)(C) ∶Aph(300 K)(Si) ∶Aph(300 K)(Ge) :Aph(300 K)(Sn) =1∶1.2∶2.2∶3.8.从图4(c)可以看出,我们所计算的非简谐性描述符Aph(300 K)与体模量BL的变化趋势一致,表明该描述符适用于同族晶体之间的非简谐性的量化比较.Ge,GaAs,CdTe,NaCl 是以Ge 为中心横向延伸的同周期元素组成的材料.我们的计算结果表明,300 K 条件下,Aph(Ge) :Aph(GaAs) ∶Aph(CdTe) ∶Aph(NaCl)=1∶2.7∶5.8∶8.7.我们所计算的晶体非简谐性描述符Aph(300 K)与晶体体模量BL变化趋势一致,同理可推论该描述符适用于同周期晶体之间的非谐性的量化比较.在计算相同元素组成但不同构型的CsPbI3时,Aph(300 K)(γ-CsPbI3) :Aph(300 K)(β-CsPbI3)=1∶2.7.目前研究认为β-CsPbI3的非简谐性高于γ-CsPbI3[57],图4(c)表明两种不同构型的CsPbI3的非简谐描述符与其材料体模量的变化趋势一致,同理可推论该描述符适用于相同元素比例下不同构型的晶体的非谐性的量化比较.

此外,非简谐性是导致材料低热导率的一个重要原因[58,7,10],最近 Knoop等[8]研究了强非简谐性在超低热导率材料中的微观作用机制,通过非简谐角度对465 种材料进行了分类,并从中找出了6 种超低热导率材料.非简谐性会限制声子的热传导,因此通常情况下晶格热导率会随着晶格非简谐性的增强而减弱.我们计算了所选9 种材料在300 K下的晶格热导率,并结合已知实验所测300 K 下的材料热导率绘制了κL(300 K) 与我们提出的非简谐描述符Aph(300 K) 之间的关系.由于计算时选取的模型是理想模型,而实际材料会存在缺陷、杂质等不可避免地影响测量数值的因素,因此理论值与实验值存在一些误差,但理论模型与实验值存在的差异并不影响材料热导率的变化趋势[59].图4(d)显示了κL(300 K)与Aph(300 K) 之间存在明显的线性关系,所选9 种材料的晶格热导率的变化趋势与我们所提出的晶体非简谐性描述符的变化趋势一致.这一线性趋势验证了我们提出的晶体非简谐性描述符Aph(T) 的准确性.

从图4(c),(d)中可以看出,两幅图中强非简谐性材料CsPbI3的非简谐性描述符偏离线性趋势的幅度均大于其余7 种弱非简谐性材料,这是因为在我们的假设中原子振动平衡位置是固定的,而类似于CsPbI3这类强非简谐性材料在高温下的原子振动平衡位置会随着振动而改变,这增大了我们所提出的描述符对于此类材料的描述误差.因此在未来的工作中将着重于探寻强非简谐性材料的原子振动平衡位置的分布规律,以改进我们所提出的材料非简谐性描述符,以求提高此描述符对强非简谐性材料的描述准确性.

4 结论

通过对周期表中以Ge 为中心及其与周围元素组成的晶体(C,Si,Ge,Sn,GaAs,CdTe,NaCl,γ-CsPbI3和β-CsPbI3)的第一性原理计算及分析,引入了一种新的描述符来量化材料的非谐性.区别于传统非简谐模拟所需要的复杂耗时的MD 计算,我们所提出的描述符只需要材料声子谱、静态自洽计算以及后续的数据拟合,大大简化了过程,降低了计算成本.将总的晶格非谐性分解成单个声子振动模式,引入声子非谐性描述符.然后,将其与玻色-爱因斯坦分布相结合,提出了温度依赖的非谐性描述符Aph(T) .计算了具有广泛代表性的9 种材料的晶体体模量和300 K 下的晶格热导率,两者与我们所提出的定量描述符Aph(300 K) 之间均具有明显线性趋势,此结果验证了所提出描述符的准确性.本研究为未来受非简谐效应驱动的特殊功能材料的高通量研究提供了范例,为其设计和应用开辟了新的可能性.此外,对于以CsPbI3为例的强非简谐材料,在高温下,其原子的平衡位置不固定,导致使用我们所提出的描述符量化非简谐性的准确性降低.为了解决这一局限性,我们未来的研究将致力于研究强非简谐材料在高温下的原子平衡位置的分布,以求提出一种更准确的理论方法来量化强非简谐材料中的非简谐性.