前言

當光通過材料時，會與材料的原子（離子）和電子相互作用，并產生反射和吸收。通過研究材料的光學性質，可以直接獲得有關電子結構、缺陷（或雜質態）以及振動模式等多方面信息。

第一性原理計算光學性質的過程就是計算介電函數的過程，本文將從表征光學性質的7個重要參數的物理意義出發，介紹它們與介電函數的關系，然后給出最簡單的密度泛函理論（DFT）加外場擾動的計算方法。最后我們會介紹PWmat計算介電函數的module并著重討論二維材料的介電函數。

介電函數與光學性質

材料的光學性質都與其介電函數ε(ω)有關，本文中的介電函數都是相對介電函數，即其相對于真空的介電函數ε0的大小。光是一種電磁波，它在介質中傳播時，相當于給介質施加了外場擾動，介質則會對擾動產生響應。考慮最簡單的情況，假設入射面在xOy平面內，光的入射方向為z方向，電場分量的偏振方向沿著x，即

介質則會響應外場變化并產生極化

相應的電位移矢量為

顯然，介電函數應該是個復數

口語中的介電函數通常指的是ε(ω)的實部，它代表相對電容率，即一個材料儲存靜電能的相對能力的大小，這個值越接近1說明材料的絕緣性越好（越接近真空）。

虛部則代表損耗因子，極化率P?隨時間的變化為

在量綱上，P?隨時間的導數代表一種電流密度j，這個電流也分為前后兩部分，極化電流和傳導電流。極化電流和電場相差一個相位π/2，電場在一個周期內做功的平均值為0，因此不消化電磁場的能量；傳導電流則具有歐姆定律的電流形式，其電導率為

單位時間內將消耗能量。因此，介電函數的虛部代表了介質消耗外場能量的能力。

我們可以通過計算得到材料的介電函數，利用介電函數求得材料的復折射率，并最終得到材料的幾個重要的光學參數：

01?復折射率（complex refractive index）

光在介質中傳播時的光速為c/n，其中n為折射率，它也應該寫成復數的形式

實部即傳統定義的折射率，虛部則代表介質內光的衰減（有時也稱虛部為“消光系數”），由折射率的定義

因此復折射率和介電函數可以相互表示

或者寫成

請注意，為了形式簡便，此處我們將常數ε0當作1處理，在實際計算中如果需要就會把這個常數補回來。

02?反射率（reflectivity）

反射光的強度Ir和入射光的強度I0之比

很顯然，它是一個0-1之間的無量綱數，常寫成百分數的形式。由于光強，在光垂直入射界面的條件下（入射角和折射角均為0）

仔細觀察可以發現，這個過程中反射系數和入射方向的正負沒有關系，即無論你是從真空到介質還是從介質到真空，反射系數都是一樣的。值得注意的是，在考慮復折射率的情況下，反射波和入射波的振幅之比也是個復數

存在輻角θ?意味著反射光和入射光之間會存在相位差，盡管這并不影響反射系數的大小，輻角θ?滿足

有時會將反射波與入射波的振幅的比值稱為反射系數

但是在黃昆先生的《固體物理學》中，反射系數（reflectance coefficient）的定義與我們這里的反射率相同，都是反射光與入射光的強度之比。這其中可能存在一些歷史原因，有待進一步考證。

03?吸收系數（absorption coefficient）

光在介質中傳播時沿著z方向的波矢

相應的，電場強度可以表示為

因此光在介質中沿著傳播方向以e指數的形式衰減，其衰減的系數為4πωk/c（多乘了系數2是因為光強正比于電場振幅的平方）。從光強的角度，經過界面反射后剩余的光在通過固體的過程中可能被吸收，具體表現為光強隨著入射深度z?衰減

其中β?就是吸收系數

在自然單位制下

需要特別注意的是，與反射系數不同，吸收系數是一個實際的物理量，它有屬于自己的量綱（通常為）。吸收系數，吸光度和吸收率互為不同的物理量！

04?透射率（transmissivity）

通過介質后光的強度It?和入射光的強度I0?的比值

它與材料的反射率R，吸收系數β和材料在入射方向的厚度d（此處為z?方向）都有關

有些時候會與透射系數（transmission coefficient）混用，但是更多的時候，透射系數的定義為透射波與入射波的振幅之比

05?吸光度（absorbance）

光線通過介質前的光強I0與通過介質后的光強It的比值的常用對數。聽起來有點繞口，這里直接給一個表達式

不難看出它其實就是透射率倒數的常用對數,可以用一個換底公式將其化簡

實際上，即使反射率高達90%，后面一項的數值也只是-2，通常我們也不會討論一個擁有90%以上的反射率的物質的吸光度。因此，常用的吸光度表達式為

06?吸收率（absorptance/absorptivity）

當光穿過一個材料的過程中，光的強度分為了反射，吸收和透射3個部分，因此材料的吸收率α，反射率R?和透射率τ?滿足以下關系

由于界面之間可能存在多重反射，實際的吸收率的形式非常復雜。僅考慮第一次入射-出射過程，則

07?發射率（emissivity）

物體表面熱輻射釋放的能量與絕對黑體熱輻射始放的能量的比值。它也是一個0-1之間的無量綱數，常寫成百分數的形式。它與溫度T?以及光的頻率ω?有關。在特定溫度T?下，總的發射率

其中uBB為黑體（blackbody, BB）的輻射能量。這本來是個熱輻射的概念，但是根據Kirchhoff法則，對于特定頻率的光，其吸收率和發射率相等，因此有

對充分厚且不透明的物體，其透射率為0，于是

?對黑體，反射率和透射率均為0，于是

對于一般情況，可以直接套用之前吸收率的計算公式?；蛘甙凑瘴墨I給出的，在上下界面都是平面且經歷過多重內部反射時

上述公式的來源可參考：

（1）10.1088/0022-3727/42/15/155412

（2）10.1088/0965-0393/22/7/075004

以上僅僅是最簡單的情況，實際的介電函數是一個二階張量，我們可以通過對稱性減少張量的矩陣元的數量，或者直接使用一個正交的晶格使整個張量只有對角元。求解材料的光學性質的關鍵在于求得復介電函數，但是這個求解過程又觸及了密度泛函理論（DFT）的一大痛點，即對外場的響應。DFT的交換關聯勢中并沒有包含動態的介電函數，這意味著第i?個態的電子激發并不會引起其它電子的狀態改變。針對上述問題，解決方案有3種：

01?直接加一個與時間有關的勢能V(r,t?)，利用劉維爾方程描述密度算符隨時間的變化。這種方法使用的波函數仍然是Kohn-Sham（KS）方程的波函數。該過程的核心在于無規相近似（random phase approximation, RPA）。

02?含時密度泛函理論（TDDFT）:含時的微擾理論，它可以更好的描述體系對外場的響應函數χ(ω,t?)。

03 GW近似：直接將動態的介電函數包含在自能項Σ?中。

在這個篇章中，我們將介紹第1種方法，在后續兩個篇章中我們將介紹TDDFT方法。GW方法過于昂貴，而且其形式已經高于DFT的框架，在未來將會有專門的篇章介紹GW以及GW+BSE。

DFT+外場擾動

在KS方程的基礎上，假設未受擾動時的哈密頓量為H0，其能量本征值和波函數滿足Bloch定理

其中l?代表能帶的指標，Ω?為使用的晶胞的體積。此處為沒有使用大家熟悉的n?作為能帶指標道個歉，因為n?在前后文使用太多次了。其密度算符ρ0滿足

其中f?代表Fermi-Dirac分布。

現在我們給該體系施加一個隨時間變化的外勢Vi，體系會產生一個響應勢Vs以屏蔽外場的擾動。此時體系的哈密頓量和密度算符會變為

其中ρ1為受到外場擾動后密度算符的改變。體系的電子氣變化滿足Liouville方程

由于擁有共同本征態的物理量的對易為0

由Poisson方程可知V和ρ1存在函數關系，為了簡化計算，可以作線性化（linearize）處理，即將二者的乘積視為高階小量，此時有

又有ρ0不隨時間變化。因此

對勢V?做Fourier變換，在第一Brillouin zone（BZ）內有

其中的矢量K為倒格矢，對化簡后的Liouville方程兩端同時求KS本征態的矩陣元有

對于右邊的第二項中有

其中的一個因子

求和之后又對r?做積分，這個過程非常繁瑣，因此引入RPA：在電子濃度足夠大時（固體中很容易滿足這個條件），電子的位置在空間的分布是雜亂無規的，在由r?決定的相位因子求和

中，若體系滿足平移對稱性，則無規則變化（q-q’≠0）的指數項之和為0，即

這也可以理解為電子濃度足夠高的情況下，相干疊加的數量級要遠大于無規相位的指數相加。剛才的表述是經典動力學中的表述，所以才有“電子濃度”這類表述，源于David在1951年關于相互作用的電子氣與等離激元的工作（Phys.Rev.82,625）中。筆者認為這種解釋比較直觀，所以在此使用。感興趣的讀者可以在李正中先生的《固體理論》的4-4中看到二次量子化表述下的量子運動方程的RPA，即將“異類”之間彼此激發產生的非線性項的求和歸零，因此RPA也可以視為一種線性化操作。在微擾理論中，RPA又被稱為環形圖（ring diagram）近似或圈圖修正，在此不繼續展開討論，想繼續進階的讀者可以看：

Phys.Rev.106,364(1957)

Phys.Rev.111,442(1958)

或者找一些多體物理的教材。經過RPA后有

設外勢Vi?隨著時間的演化為

假設它對體系的影響是絕熱的

因此ρ0和Vs都有充足的時間對Vi做出響應，它們對時間的變化形式均與Vi相同

于是Liouville方程的左端可以變成

最終我們得到密度ρ1在KS本征態下的矩陣元

其中

體系對V?的響應表現在極化過程，極化會帶來電荷密度的改變（束縛電荷），根據介質中的Maxwell方程組，電荷分為束縛電荷ρb和自由電荷ρf

它們可以分別用與極化矢量P?和和電位移矢量D?的散度來表示

于是，在Gaussian單位制下

從另一個角度理解上述關系，電荷密度的變化完全來源于為了屏蔽外勢Vi而產生的極化。因此，體系響應外場而產生的屏蔽勢Vs和電荷密度的變化可以有一個屬于自己的Poisson方程

其中電荷密度的變化恰好與ρ1在KS本征態下的矩陣元有關

它與V通過ρ1的矩陣元而產生直接的數學關系。于是，我們可以在Vs和V之間建立一種自洽關系

其中

這也是一個矩陣，且矩陣元與倒格矢K，K’?有關。在單電子近似下，極化率擁有如下形式

其中是k點歸一化后權重，它們的和為1。于是我們可以看出，矩陣具有以下意義

相應的，介電函數的也可以寫成與倒格矢K，K’?有關的矩陣元

為了求解介電函數的矩陣，仍然需要繁瑣的迭代，人們通常只考慮長波極限下的情況。長波極限對應的頻率較高，比較符合可見光波段的吸收過程。當|q|遠小于|K|時，可以近似的認為標量場f?在一個晶胞范圍內的平均值就是標量場在K=0時的值，且f?的平均值的梯度就是f?梯度的平均值，于是

又因為

可得

最終我們得到了介電函數的倒數

這樣得到的介電函數稱為“宏觀介電函數”

假如介電函數的矩陣沒有非對角元，則可以直接求逆

由復數的展開式

以及二階張量的表現形式

宏觀介電函數的虛部可以寫作

在只有垂直躍遷的情況下，l’?即為導帶CB，l?即為價帶VB，再引入原子單位制，則虛部可以改寫為

實部則可以直接由虛部通過Kramers-Kronig(KK)變換得到

介電函數矩陣經常會存在不為0的非對角元，這是由于微觀上單個偶極子處的電場強度并不是平均的宏觀電場的強度，而應該是在單個電偶極子中心處的平均電場強度Eeff，這是一種局域場強。Eeff和宏觀電場的差別來源正是介質中電偶極子的相互作用。在晶體的計算中，這種局域場效應則來源于周期性的勢場。有多種方法可以修正局域場效應，我們將在今后的極化性質篇章，結合離子部分的介電函數等因素，詳細討論這個問題。

本部分的推導用到但是不限于以下參考文獻：

1.???? Phys.Rev.115, 786(1959)

2.???? Phys.Rev.129, 62(1963)

3.???? Phys.Rev.B.33,7017(1986)

PWmat中利用RPA計算光學性質的module

PWmat具有兩個module可以利用RPA計算光學性質：module 18的前一部分和module 38。它們的前置工作是類似：

01?做自洽計算，得到波函數。

可以選擇在這一步就將k-mesh取得很密，同時保證足夠的空帶數目。也可以選擇在這一步使用較少的k-mesh點和空帶數，然后在自洽計算的基礎上做一個加密K點和空帶數目的非自洽計算并輸出波函數。

02?利用第一步的波函數做DOS計算。

兩個module的區別在于對DOS的插值方式不同。計算光學性質需要極密的k-mesh，PWmat提供了兩種插值方式，第一種是在每個k點周圍以立方體的形式插值，第二種則是在每個k點周圍以四面體的形式插值。計算種，推薦使用第二種插值方式。與module 18相比，module 38不僅考慮了局域場效應，還增加了贗勢中非局域部分帶來的影響。

在實際計算中，可以通過平移算符把對k?的小位移移出Bloch波函數，這樣就可以用普通的KS方程的本征態來表示介電函數的虛部

上述展開用了長波極限下的一個近似

由于有周期性邊界條件的存在，在晶體中直接定義位移r?的矩陣元是有問題的。解決方法是將它與晶格的自洽場哈密頓量HSCF一起考慮

勢能中局域的那一部分剛好與位移r?對易，因此哈密頓量只會剩下動能與r?的對易子和非局域勢Vnl與r?的對易子

以此修正非局域勢對介電函數的影響。

PWmat中有JOB = MOMENT可以進行上述計算，在算完DOS后再計算MOMENT, 然后使用腳本RPA_absorb.x即可獲得介電函數。在20220401的版本更新后，這個腳本還會輸出復折射率，反射率和消光系數。用戶也可以通過腳本RPA_diel_G1G2_000.x來計算宏觀靜態介電函數。下圖展示的是集成module 38的GUI的計算結果（龍訊超算云平臺Mcloud介紹）。

想要計算發射譜的用戶則可以參照module 55。

彩蛋：二維材料的介電函數

在計算二維材料的性質時，為了能使用Bloch波函數，我們會在垂直方向添加真空層，構成一個超胞。于是當我們計算二維材料的介電函數時，實際上是把真空層也當成了bulk的一部分進行計算，這樣直接得到的結果肯定是錯誤的?？紤]到二維材料在面內仍然是正常的晶體，由介電函數的虛部和實部的形式

對面內來說，介電函數的錯誤來自于晶胞的體積取大Ω了，于是面內的函數算出來是偏小的。假設二維材料的厚度為t，真空層方向的晶格常數為c，面內的介電函數與計算所得的介電函數的關系是

至于如何定義二維材料的厚度，則見仁見智?？梢灾苯佑米罡吆妥畹偷脑觶?坐標之差來定義，也可以對z方向的local potential做個積分，去掉能量為真空能級的部分，剩下的部分即為層厚。如果在寫paper的時候想偷懶，則做到這一步就可以結束了。你可以說：“假設我的入射光方向是沿著z方向，入射光的偏振方向都在xOy平面內，所以我要計算的所有光學性質只需要面內的分量?！彼追Q語言的藝術（耍流氓）。

假如遇到一個認真的審稿人，偏要讓入射光斜著打；又或者你本身有比較高的學術追求，很想弄清楚垂直層方向的介電函數究竟是什么樣的。我們在此提供一個較為簡單的思路，參考文獻的DOI如下：

1.???? 10.1039/c8cp04743j

2.???? 10.1038/s41699-018-0050-x

3.???? 10.1038/s41467-021-25310-2

這種思路的出發點是把垂直層的方向視為兩個電容的串聯：二維材料和真空。根據串聯電容的表達式

于是可以得到

這是宏觀靜態介電函數的修正方式，它也可以應用于介電函數的實部，而虛部可以近似使用面內的修正方式。

更進一步，我們可以嘗試將實部和虛部直接引入電容的關系中，解如下的方程組

歡迎感興趣的讀者測試并反饋結果。

---

公司簡介
北京龍訊曠騰科技有限公司是成立于2015年4月的國家高新技術企業，是國內材料計算模擬工具軟件研發創新的領導者，致力于開發滿足“工業4.0”所需的原子精度材料研發Q-CAD（quantum-computer aided design）軟件。公司自主開發的量子材料計算軟件PWmat（平面波贗勢方法并基于GPU加速）可以進行電子結構計算和從頭算分子動力學模擬，適用于晶體、缺陷體系、半導體體系、金屬體系、納米體系、量子點、團簇和分子體系等。官網為www.pwmat.com，微信公眾號為：PWmat1。

總結

以上是生活随笔為你收集整理的【龙讯module小课堂】“光”怪陆离：PWmat计算光学性质（一）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。