材料的非簡諧性描述符*

李環婭 周柯 尹萬健

(蘇州大學能源學院,能源與材料創新研究院,江蘇省先進碳材料與可穿戴能源技術重點實驗室,蘇州 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為例的強非簡諧材料,在高溫下,其原子的平衡位置不固定,導致使用我們所提出的描述符量化非簡諧性的準確性降低.為了解決這一局限性,我們未來的研究將致力于研究強非簡諧材料在高溫下的原子平衡位置的分布,以求提出一種更準確的理論方法來量化強非簡諧材料中的非簡諧性.