# Fileset

[原稿.pdf](https://mdr.nims.go.jp/filesets/ef7001dc-8395-430c-a3ef-74482dc11e1b/download)

## Creator

[只野 央将](https://orcid.org/0000-0002-8132-2161)

## Rights

[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[第一原理計算を用いた熱伝導・電気伝導予測の前線](https://mdr.nims.go.jp/datasets/7e8552e7-ab2e-4e4a-bfdb-90eec461f732)

## Fulltext

1．はじめに   量子論に基づく電子状態計算，特に密度汎関数理論（DFT）に基づくバンド計算はセラミックス材料を含むあらゆる材料物性を解析するツールとして広く普及している．近年では大量の DFT計算の結果から機械学習モデルを構築し，それによって効率的に新材料設計を行う試み 1)が報告されており，予測の信頼性には改善余地があるものの 2)，電子状態計算の存在感は以前にも増している．DFT計算を行えばバンド構造や電子状態など基本的な物性は比較的容易に解析可能である．一方，実際の材料設計で興味がある複素誘電率などの光学特性や，超伝導転移温度，電気・熱伝導率などを非経験的に予測するのはより難しく，これまでは現象論に基づく解析に頼らざるを得なかった．しかし，効率的な計算手法の確立やアクセス可能な計算資源の増加のおかげで，現代ではこれらの物性値を第一原理的に定量測出来るようになっている．本稿では、熱電変換材料の設計で重要な熱伝導率，電気伝導率，およびゼーベック係数の第一原理計算に焦点を当て，近年の計算手法の進展と逆ペロブスカイト Ba3SiO, Ba3GeOへの応用例 3)を紹介し，その有効性を示す． 2．熱伝導率   固体における熱伝導の担い手は格子振動（フォノン）と荷電粒子であり，磁性体の場合はスピン励起（マグノン）による寄与もある．セラミックス材料や熱電材料はバンドギャップが開いているか，あるいは金属的伝導を示す場合でもキャリア濃度がそれほど高くない場合が多く，そこではフォノンによる熱伝導が支配的になる．フォノン熱伝導を第一原理から評価する場合，Boltzmann 方程式が広く用いられており，緩和時間近似（RTA）の元で以下の式から計算される． 𝜅ph = 1𝑁𝑞 ∑ 𝑐𝒒ν𝒒ν𝒗𝒒𝜈 ⊗ 𝒗𝒒𝜈τ𝒒ν (1) ここで𝒒, νはフォノンの運動量ベクトルとバンドインデックスであり，𝑐𝒒𝜈 , 𝒗𝒒ν, τ𝒒νはそれぞれフォノンモード𝒒νの比熱，群速度，および緩和時間である．比熱と群速度はそれぞれフォノン分散曲線から得られる一体の物理量であるが，緩和時間τ𝒒νを決定するにはフォノン散乱効果を取り入れる必要がある．フォノン散乱過程には，格子振動の非調和性に起因するフォノン−フォノン散乱，電子格子相互作用に起因するフォノン−電子散乱，さらにフォノン−不純物散乱やフォノン−粒界散乱がある．このうち前者二項は内因性の散乱過程であり，第一原理的に評価可能である．一方，外因性の散乱過程については汎用的な評価手法が確立していないが，同位体不純物によるフォノン散乱確率に関しては古くに提案された計算手法 4)が利用されており，また異種元素置換や空孔によるフォノン散乱についてもいくつか計算手法が提案されている 5),6),7)．これらの散乱効果をすべて考慮できれば理想的だが，実際には低温域を除く温度域で主要項となるフォノン−フォノン散乱を考慮すれば十分であることが多い．その中で最も重要なのは，3 つのフォノンが関わる 3 フォノン散乱過程であり，その散乱強度は Π𝒒𝜈(3)(𝜔) = 116𝑁𝑞 ∑ℏ|Φ(3)(−𝒒𝜈; 𝒒1𝜈1; 𝒒2𝜈2)|2𝜔𝒒𝜈𝜔𝒒1𝜈1𝜔𝒒2𝜈2𝒒1𝜈1,𝒒2𝜈2  ×  Δ(−𝒒 + 𝒒1 + 𝒒2)𝑓 (1, 2, 𝜔𝑐) (2) の虚部で与えられる．Π𝒒𝜈(3)(𝜔)はフォノンの自己エネルギーと呼ばれる量である．ここで，Φ(3)(−𝑞𝜈; 𝑞1𝜈1; 𝑞2𝜈2)は 3 つのフォノンモードの相互作用係数であり，これはポテンシャルエネルギーの 3 階微分である 3 次非調和力定数（IFC）から計算する．Δ(𝒒)は𝒒が逆格子ベクトルの整数倍の時に 1，それ以外の時にゼロとなる関数であり，運動量保存を反映している．一方，𝑓 (1, 2, 𝜔𝑐)はエネルギー保存を反映した項であり， 𝑓 (1,2, 𝜔𝑐) = ∑ 𝜎 [1 + 𝑛1 + 𝑛2𝜔𝑐 + 𝜎(𝜔1 + 𝜔2)−𝑛1 − 𝑛2𝜔𝑐 + 𝜎(𝜔1 − 𝜔2)] 𝜎=±1 で定義される．ここで 𝑛𝑖 = [exp (ℏ𝜔𝑖/𝑘𝐵𝑇 ) − 1]−1 はBose-Einstein 分布関数，そして𝜔𝑐 = 𝜔 + 𝑖0+ (0+は正の微少量)である．式(2)から求まるフォノン散乱確率𝜏𝒒𝜈−1 = (2/ℏ)ImΠ𝒒𝜈(3)(𝜔𝒒𝜈)は Fermi の黄金律による導出と等価である．式(2)で温度に依存する項は𝑓 (1, 2, 𝜔𝑐)に含第一原理計算を用いた熱伝導・電気伝導予測の前線  Frontier of thermal and electrical conductivity prediction using first-principles calculations Keywords: First-principles calculation, phonons, electrons, thermoelectricity, inverse perovskite 只野 央将 Terumasa TADANO (National Institute for Materials Science) まれる𝑛𝑖項のみで，これは高温極限で温度𝑇の線形関数になる．したがって，高温域で𝜏𝒒𝜈−1 ∝ 𝑇となり，式(1)の比熱𝑐𝒒𝜈は高温域でほぼ Dulong-Petit 極限で一定になることを考えると𝜅ph ∝ 𝑇 −1が得られる．この温度依存性は多くの単結晶材料で見られるものと整合する．以上の式(1)，(2)にもとづく熱伝導計算が最も標準的な第一原理アプローチであるので，以下では「従来法」と呼ぶことにする．従来法は，例えば ShengBTE, phono3py, ALAMODE, Phoebe などさまざまなソフトウェアで利用可能であり，共有結合性が強い半導体や絶縁体を中心に良好な予測精度を示す事が確認されている 8)．では以下では従来法が不十分なケースとその場合に用いられるより先進的な計算手法を紹介しよう． 一つ目は，4 つのフォノンが関わる 4 フォノン散乱の効果で，この重要性が初期に示されたのはシリコンや BAsなどの高熱伝導材料である．特に BAsにおいてはその影響は大きい．BAs は，従来法で予測された極めて高い室温熱伝導（~2,300 W/mK，ダイヤモンドに匹敵）が当初注目を浴びたが 9)，その後の単結晶での熱伝導測定の結果，熱伝導率は予想の半分程度(~1,200 W/mK)しか無かった 10)．この差を説明するには，式(2)に加えて 4フォノン散乱過程を考慮し，フォノン散乱確率の評価式を𝜏𝒒𝜈−1 = (2/ℏ)[ImΠ𝒒𝜈(3)(𝜔𝒒𝜈) + ImΠ𝒒𝜈(4)(𝜔𝒒𝜈)]と改める必要があると結論されている．Π𝒒𝜈(4)(𝜔𝒒𝜈)の具体的な計算式は他文献 11)に譲るが，Π𝒒𝜈(4)(𝜔𝒒𝜈)の計算には 4 次の非調和 IFCが必要となり，また積分自体もさらに多重になるため計算コストが大きく跳ね上がる．なお，Π𝒒𝜈(4)(𝜔𝒒𝜈)は高温極限で𝑇 2の温度依存性を示すため，4 フォノン散乱を考慮する場合の熱伝導率は𝜅ph ∝𝑇 −𝛼 (𝛼 > 1)となり従来法よりも温度増加に伴う低下が速くなる．しかし，熱伝導率が低い熱電材料などでは，熱伝導率の温度依存性が𝜅ph ∝ 𝑇 −1よりも緩やかな事が多いため，4 フォノン散乱を考慮すれば実験結果をより良く再現するとも限らない点に注意が必要である． 二つ目はフォノン振動数や分極ベクトルの有限温度効果で，これはペロブスカイトなど構造相転移を起こす物質系では特に重要になる．例えば 105 Kで正方晶−立方晶転移を示す SrTiO3では，構造相転移に関連する低エネルギーモード（ソフトモード）の振動数が顕著な温度依存性を示すが，従来法は原子変位が格子間隔に比べて十分小さい極限で正当化される調和近似を用いており，フォノンの温度変化が記述できない．これを解決するには，有限の熱振幅下で非調和効果を取り込めるフォノン計算法を用いる必要がある．それを実現する方法の一つが自己無撞着フォノン（SCP）法 12)である．SCP法では原子座標の熱分布が Gaussianであると仮定し，その元で非調和効果を平均場近似で取り入れる．より具体的には，以下の式を解けば良い． Φ𝑖𝑗eff (𝑇 ) = Φ𝑖𝑗 + 12 ∑ Φ𝑖𝑗𝑘𝑙𝑘𝑙〈𝑢𝑘𝑢𝑙〉eff  . (3) ここで，Φ𝑖𝑗およびΦ𝑖𝑗𝑘𝑙は 2 次および 4 次 IFCであり，𝑢𝑖は平衡点からの原子変位量である．左辺のΦ𝑖𝑗eff (𝑇 )は非調和効果によって補正された有効的な 2 次 IFCであり，これを用いて動力学行列を構成し，それを対角化することで有限温度でのフォノン振動数Ω𝒒𝜈(𝑇 )が求まる．一方，右辺の〈𝑢𝑘𝑢𝑙〉effは温度に依存する平均自乗変位である．このアンサンブル平均は左辺Φ𝑖𝑗eff (𝑇 )の結果から計算するため，式（3）は両辺を自己無撞着に解く必要がある．筆者らは SCP法を立方晶 SrTiO3へ適用することで，従来法では不可能だった高温相のフォノン計算および熱伝導計算を実現した 12)．そこでは式(3)を解くことで，さまざまな温度でのフォノン振動数Ω𝒒𝜈(𝑇 )および分極ベクトルを計算し，それを式(1), (2)の𝜔𝒒𝜈のかわりに使うことでフォノン散乱強度と𝜅phを求めた．この手続きによって，𝜅ph ∝ 𝑇 −𝛽(𝛽~0.6 − 0.7)という𝜅ph ∝ 𝑇 −1よりも緩やかな温度依存性まで含めて実験結果を定量的に説明する事に成功した．温度依存性が𝜅ph ∝ 𝑇 −1よりも緩やかである理由は以下のように説明できる．非調和効果によってフォノン振動数は温度上昇とともにハード化（増加）する．この傾向は低エネルギー領域の光学フォノンで特に顕著である．すると，𝑓 (1, 2, 𝜔𝑐)に含まれる𝑛𝑖の温度変化が𝑛𝑖 ∝ 𝑇よりも緩やかになり，さらに光学フォノンと音響フォノンの混成も弱まるため，式（2）の散乱強度が𝜏𝒒𝜈−1 =(2/ℏ)ImΠ𝒒𝜈(3)(𝜔𝒒𝜈) ∝ 𝑇よりも穏やかになる．したがって，式(1)の熱伝導率も𝜅ph ∝ 𝑇 −𝛽(𝛽 < 1)となる．SCP法は熱伝導シミュレーションだけでなく，有限温度での変位型構造相転移 13)や熱膨張 14)の効率的なシミュレーションにも応用可能であり，SCP法を利用した研究例は増加傾向にある．また，SCP 法を超える試み 15),16)もなされており，今後も更なる理論的な進展が期待される．  以上の他にも，フォノン散乱が極めて強く，フォノン準粒子描像など Boltzmann 輸送方程式の前提が破綻する場合にも有効な計算手法が提案されるなど 17),18)，多くの基礎的な進展が見られる．これらの効果をすべて考慮すれば理論的には最も高い予測精度が期待されるが，そのレベルの計算が常に実験結果を定量的に説明するとも限らないところが難しい点である．例えば，従来法で無視する複数の効果が違いにキャンセルして，たまたま実験の熱伝導率を良く再現することもある．また，すべての計算の元になる DFTにもさまざまなレベルの近似が存在し，フォノン振動数や熱伝導率の値もそれに応じて変化する．したがって，定量性については第一原理計算に過度な期待はできない．しかし，従来法よりも信頼性が高い予測が可能になった点に加え，熱伝導の物質依存性などを微視的に理解する上でこれらの第一原理手法が有用である点に異論は無い． 3．電気伝導率・ゼーベック係数 熱電変換材料で熱伝導率と同じく重要になるのは電気伝導率𝜎およびゼーベック係数𝑆である．これらもBoltzmann方程式を用いて計算するのが一般的であり，𝜎 = 𝑒2𝐾0および𝑆 = − 1𝑒𝑇 𝐾0−1𝐾1から評価する．ここで，𝐾𝑗は以下で定義される． 𝐾𝑗 = 1𝑁𝑘 ∑ 𝜏𝒌𝑛𝒗𝒌𝑛𝒌𝑛⊗ 𝒗𝒌𝑛(𝜀𝒌𝑛 − 𝜇)𝑗(− 𝜕𝑓𝜕𝜀) . (4) ここで，𝜀𝒌𝑛, 𝒗𝒌𝑛, 𝜏𝒌𝑛はそれぞれ固有状態𝒌𝑛の一体バンドエネルギー，群速度，および緩和時間である．また，𝜇は化学ポテンシャル，𝑓𝒌𝑛 = {exp [(𝜀𝒌𝑛 − 𝜇)/𝑘𝐵𝑇 ] +1}−1は Fermi-Dirac 分布関数である．従来の第一原理計算では電子の緩和時間𝜏𝒌𝑛を定量計算せず，𝜏𝒌𝑛 ≈ 𝜏と近似する一定緩和時間近似（CRTA）がよく用いられてきた．例えば BoltzTraP コードでは CRTA が採用されている．この近似を用いれば式(4)の𝜏𝒌𝑛が和の外側へ移動でき，さらに𝑆では分母と分子で𝜏がキャンセルするので𝜏に依らない評価が可能である．しかし，𝜎ではそのようなキャンセルは発生しないので，定量予測には𝜏の情報が必要となる．過去の研究では𝜏の値として𝜏 ∼10 fs が近似的に用いられてきたが，それが広い材料系で良い近似になる保証は全くない． 近年では，電子格子相互作用を考慮して𝜏𝒌𝑛を第一原理的に評価することが可能になっている．具体的には電子自己エネルギーの Fan-Migdal 項 19) Σ𝒌𝑛FM(𝜔) = ∑ |𝑔𝑚𝑛𝜈(𝒌, 𝒒)|2𝑚,𝒒𝜈          × [𝑛𝒒𝜈 + 𝑓𝒌+𝒒𝑚𝜔 − 𝜀𝒌+𝒒𝑚 + 𝜔𝒒𝜈 − 𝑖0+ + 𝑛𝒒𝜈 + 1 − 𝑓𝒌+𝒒𝑚𝜔 − 𝜀𝒌+𝒒𝑚 − 𝜔𝒒𝜈 − 𝑖0+] , (5) の虚部から𝜏𝒌𝑛−1 = (2/ℏ)Im Σ𝒌𝑛FM(𝜀𝒌𝑛)と求める．ここで，𝑔𝑚𝑛𝜈(𝒌, 𝒒) = ⟨𝜓𝒌+𝒒𝑚|𝜕𝒒𝜈𝑉scf |𝜓𝒌𝑛⟩は電子格子結合定数であり，原子座標の微小摂動に対する DFT計算で用いる有効ポテンシャル𝑉scfの応答から計算される量である．これは密度汎関数摂動論（DFPT）を用いて効率的に計算可能である． 電気伝導率やゼーベック係数には，式(4)からも分かるとおり，フェルミエネルギー𝜀𝐹（または化学ポテンシャル）から測って±𝑘𝐵𝑇程度の狭いエネルギー領域の状態のみ寄与する．したがって，計算結果を収束させるために非常に密に𝒌点サンプリングをする必要がある．しかし，密な𝒌点で𝑔𝑚𝑛𝜈(𝒌, 𝒒)を DFPT計算するのはコスト的に困難である．これに対処するため，最局在Wannier関数を用いた補間方法が広く利用されている．Wannier 補間を用いれば，比較的疎な𝒌点と𝒒点で𝑔𝑚𝑛𝜈(𝒌, 𝒒)を DFPT計算し，それを何倍も密な𝒌, 𝒒点へ精密に補間出来る 19)．この際，輸送に寄与しうる𝜀𝐹近傍のバンドを正確に再現する最少の Wannier 軌道の組を利用すれば，補間のコストも大幅に削減できる．しかし，最適な最局在 Wannier 関数を構築するのは多くの場合簡単ではなく，ある程度の試行錯誤が必要である． 近年では，上述の DFPT+Wannierよりも簡便な𝜏𝒌𝑛の計算手法も利用可能である．AMSET コード 20)では電子散乱過程として，音響フォノンによる変形ポテンシャル散乱，イオン化不純物散乱，極性光学フォノン散乱，そしてピエゾ散乱を考え，それぞれの現象論に含まれるパラメータを第一原理計算から決定するというアプローチを採用している．この方法は，さまざまな近似を用いているため精度は DFPT+Wannierに劣るが，計算コストが非常に低いのが魅力である．さらに，広い物質系に対して DFPT+Wannierと定量的にも近い輸送係数を与える事が報告されており 20)，CRTA を超える第一歩として魅力的な選択になっている．ただし，実験結果に頼らずに新物質の𝜎, 𝑆や熱電変換性能𝑍𝑇を理論予測する場合には，利用する近似が対象物質において妥当なものか十分吟味することを推奨する． 4．逆ペロブスカイト物質への応用とまとめ これまで解説してきた先進手法を Ba3BO (B=Si, Ge)へ適用した最近の研究例を紹介し，手法の有効性を議論する．Ba3BOは比較的高い p 型の熱電変換性能𝑍𝑇を有することが近年実験的に見出された 3)．結晶構造はいわゆる逆ペロブスカイト構造であり，酸素原子の周りを Ba原子が八面体配位している（図 1a）．室温で空間群は Pnmaである．多結晶試料で熱伝導率を測定したところ，室温で𝜅 ∼ 0.77 − 1.00 W/mKという極めて低い値が得られた．これは構成元素が似たペロブスカイト SrTiO3と比較して一桁小さい．また，パワーファクター𝜎𝑆2も十分大きな値が実現され，結果として温度 600 Kで 0.8を超える𝑍𝑇値が Ba3SiOで達成された．そこで，この低熱伝導の起源が何か，さらに𝑍𝑇値を向上する余地はあるのかという点を解明するため，第一原理計算を行った．フォノン物性にはALAMODE コードを利用し，SCP計算と 3フォノン散乱を考慮した熱伝導率計算を実施した．また，𝜏𝒌𝑛の計算には DFPT+Wannier法を利用し，電子輸送係数は Perturbo コード 21)により求めた．その際，孤立した12本の valenceバンド（図 1b）を再現するために 12個の𝑝軌道を用いて最局在Wannier関数を構築した．計算条件の更なる詳細については原論文 3)を参照されたい． 室温でのフォノン分散曲線を図 1cに示す．SrTiO3と比較すると Ba3BOのフォノン振動数は大きく低下しているのが見える．これは八面体を構成する Ti–O結合に比べて O–Ba結合の共有結合性が弱いためと理解できる．その結果，音響フォノン群速度も SrTiO3の半分程度となっている（図 1d）．また，SrTiO3で顕著な群速度が大きな光学フォノンモードも Ba3BOには存在しない．一方，同じエネルギーで比較すれば，フォノン緩和時間τ𝒒νは SrTiO3と大きな違いは無い（図 1d）．したがって，低熱伝導率の主な原因は，共有結合性の低い O–Ba ネットワークによるフォノン群速度の低下であると結論づけられる．  図 2に熱伝導率および𝑍𝑇値の実験結果と計算結果の比較を示す．第一原理計算によって低熱伝導を定量的に再現できている（図 2a）．また，𝑍𝑇値についても，制御可能な計算パラメータが存在しないことを考えれば，実験結果を十分に再現できていると言える（図 2b）．なお，Ba3GeOで 500 K以上で見られる𝑍𝑇値の減少は，構造相転移の影響があるものと推察されている 3)．最後に，第一原理計算から予測した𝑍𝑇値のキャリア濃度依存性を図 2cに示す．この予測によると，キャリア濃度を 1020 cm-3程度まで増やすことが出来れば Ba3SiOの𝑍𝑇値がさらに上昇し，600 Kで𝑍𝑇＝2に到達し得ることを示している．この予測の実験的な検証は今後の課題だが，本成果は環境負荷の低い元素のみから構成される逆ペロブスカイトが熱電材料としての高いポテンシャルを持つことを示しており，組成最適化など今後のさらなる進展が期待される．  本稿では第一原理計算を用いた熱伝導・電気伝導およびゼーベック係数の予測手法について，近年の発展を中心に解説を行った．また，その手法がどのように物性の理解や材料設計に応用できるかを，Ba3BOを例に実演した．本稿が熱・電子輸送特性の第一原理計算を始めるきっかけになれば幸いである．   謝  辞 本稿で紹介した研究の一部は片瀬貴義，HE Xinyi，神谷利夫の各氏をはじめとする共同研究者との研究成果に基づきます．この場を借りてお礼申し上げます．     図 1 Ba3BO の結晶構造と電子バンド構造(Ba3SiO)，および・フォノン分散曲線・群速度・緩和時間の SrTiO3との比較．文献 3)より CC BY 4.0 のもとレイアウトを一部改変し掲載．  a bc d 図 2 Ba3BO の熱伝導率，𝑍𝑇値の温度依存性，および計算から予測された𝑍𝑇値のキャリア濃度依存性．文献 3)より CC BY4.0のもとレイアウトを一部改変し掲載．  a bc文  献 1)  A. Merchant, S. Batzner, S. S. Schoenholz, M. Aykol, G. Cheon, and E. D. Cubuk, Nature 624, 80 (2023). 2)  A. K. Cheetham and R. Seshadri, Chem. Mater. 36, 3490 (2024).  3)  X. He, S. Kimura, T. Katase, T. Tadano, S. Matsuishi, M. Minohara, H. Hiramatsu, H. Kumigashira, H. Hosono, and T. Kamiya, Adv. Sci. e2307058 (2023). 4) S.-I. Tamura, Phys. Rev. B 27, 858 (1983). 5) N. A. Katcho, J. Carrete, W. Li, and N. Mingo, Phys. Rev. B 90, 094117 (2014). 6)  F. Körmann, Y. Ikeda, B. Grabowski, and M. H. F. Sluiter, npj Comput. Mater. 3, 36 (2017). 7) S. Thébaud, C. A. Polanco, L. Lindsay, and T. Berlijn, Phys. Rev. B 102, 094206 (2020). 8)  L. Lindsay, Nanoscale Microscale Thermophys. Eng. 20, 67 (2016).  9)  L. Lindsay, D. A. Broido, and T. L. Reinecke, Phys. Rev. Lett. 111, 025901 (2013). 10) Tian, B. Song, X. Chen, N. K. Ravichandran, Y. Lv, K. Chen, S. Sullivan, J. Kim, Y. Zhou, T.-H. Liu, M. Goni, Z. Ding, J. Sun, G. A. G. Udalamatta Gamage, H. Sun, H. Ziyaee, S. Huyan, L. Deng, J. Zhou, A. J. Schmidt, S. Chen, C.-W. Chu, P. Y. Huang, D. Broido, L. Shi, G. Chen, and Z. Ren, Science 361, 582 (2018). 11) Z. Han, X. Yang, W. Li, T. Feng, and X. Ruan, Comput. Phys. Commun. 270, 108179 (2022). 12) T. Tadano and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015). 13) R. Masuki, T. Nomoto, R. Arita, and T. Tadano, Phys. Rev. B 106, 224104 (2022). 14) R. Masuki, T. Nomoto, R. Arita, and T. Tadano, Phys. Rev. B. 107, 134119 (2023).  15) T. Tadano and W. A. Saidi, Phys. Rev. Lett. 129, 185901 (2022). 16) E. Xiao and C. A. Marianetti, Phys. Rev. B 107, 094303 (2023).  17) M. Simoncelli, N. Marzari, and F. Mauri, Nat. Phys. 15, 809 (2019).  18) L. Isaeva, G. Barbalinardo, D. Donadio, and S. Baroni, Nat. Commun. 10, 3853 (2019).  19) F. Giustino, Rev. Mod. Phys. 89, 015003 (2016).  20) A. M. Ganose, J. Park, A. Faghaninia, R. Woods-Robinson, K. A. Persson, and A. Jain, Nat. Commun. 12, 2222 (2021). 21) J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M. Bernardi, Comput. Phys. Commun. 264, 107970 (2021).