米乐m6官方网站:【龙讯module小课堂】“光”怪陆离：PWmat计算光学性质（一）
当光通过材料时，会与材料的原子（离子）和电子相互作用，并产生反射和吸收。通过研究材料的光学性质，可以直接获得有关电子结构、缺陷（或杂质态）以及振动模式等多方面信息。
第一性原理计算光学性质的过程就是计算介电函数的过程，本文将从表征光学性质的7个重要参数的物理意义出发，介绍它们与介电函数的关系，然后给出最简单的密度泛函理论（DFT）加外场扰动的计算方法。最后我们会介绍PWmat计算介电函数的module并着重讨论二维材料的介电函数。
材料的光学性质都与其介电函数 有关，本文中的介电函数都是相对介电函数，即其相对于真空的介电函数 的大小。光是一种电磁波，它在介质中传播时，相当于给介质施加了外场扰动，介质则会对扰动产生响应。考虑最简单的情况，假设入射面在xOy平面内，光的入射方向为z方向，电场分量的偏振方向沿着x，即
口语中的介电函数通常指的是 的实部，它代表相对电容率，即一个材料储存静电能的相对能力的大小，这个值越接近1说明材料的绝缘性越好（越接近真空）。
在量纲上，P 随时间的导数代表一种电流密度j，这个电流也分为前后两部分，极化电流和传导电流。极化电流和电场相差一个相位π/2，电场在一个周期内做功的平均值为0，因此不消化电磁场的能量；传导电流则具有欧姆定律的电流形式，其电导率为
我们可以通过计算得到材料的介电函数，利用介电函数求得材料的复折射率，并最终得到材料的几个重要的光学参数：
实部即传统定义的折射率，虚部则代表介质内光的衰减（有时也称虚部为“消光系数”），由折射率的定义
请注意，为了形式简便，此处我们将常数ε0当作1处理，在实际计算中如果需要就会把这个常数补回来。
很显然，它是一个0-1之间的无量纲数，常写成百分数的形式。由于光强I∝Ex2，在光垂直入射界面的条件下（入射角和折射角均为0）
仔细观察可以发现，这个过程中反射系数和入射方向的正负没有关系，即无论你是从真空到介质还是从介质到真空，反射系数都是一样的。值得注意的是，在考虑复折射率的情况下，反射波和入射波的振幅之比也是个复数
存在辐角θ 意味着反射光和入射光之间会存在相位差，尽管这并不影响反射系数的大小，辐角θ 满足
但是在黄昆先生的《固体物理学》中，反射系数（reflectance coefficient）的定义与我们这里的反射率相同，都是反射光与入射光的强度之比。这其中可能存在一些历史原因，有待进一步考证。
因此光在介质中沿着传播方向以e指数的形式衰减，其衰减的系数为4πωk/c（多乘了系数2是因为光强正比于电场振幅的平方）。从光强的角度，经过界面反射后剩余的光在通过固体的过程中可能被吸收，具体表现为光强随着入射深度z 衰减
需要特别注意的是，与反射系数不同，吸收系数是一个实际的物理量，它有属于自己的量纲（通常为cm-1）。吸收系数，吸光度和吸收率互为不同的物理量！
有些时候会与透射系数（transmission coefficient）混用，但是更多的时候，透射系数的定义为透射波与入射波的振幅之比
实际上，即使反射率高达90%，后面一项的数值也只是-2，通常我们也不会讨论一个拥有90%以上的反射率的物质的吸光度。因此，常用的吸光度表达式为
当光穿过一个材料的过程中，光的强度分为了反射，吸收和透射3个部分，因此材料的吸收率α，反射率R 和透射率τ 满足以下关系
物体表面热辐射释放的能量与绝对黑体热辐射始放的能量的比值。它也是一个0-1之间的无量纲数，常写成百分数的形式。它与温度T 以及光的频率ω 有关。在特定温度T 下，总的发射率
其中uBB为黑体（blackbody, BB）的辐射能量。这本来是个热辐射的概念，但是根据Kirchhoff法则，对于特定频率的光，其吸收率和发射率相等，因此有
对于一般情况，可以直接套用之前吸收率的计算公式。或者按照文献给出的，在上下界面都是平面且经历过多重内部反射时
以上仅仅是最简单的情况，实际的介电函数是一个二阶张量，我们可以通过对称性减少张量的矩阵元的数量，或者直接使用一个正交的晶格使整个张量只有对角元。求解材料的光学性质的关键在于求得复介电函数，但是这个求解过程又触及了密度泛函理论（DFT）的一大痛点，即对外场的响应。DFT的交换关联势中并没有包含动态的介电函数，这意味着第i 个态的电子激发并不会引起其它电子的状态改变。针对上述问题，解决方案有3种：
①直接加一个与时间有关的势能V(r,t )，利用刘维尔方程描述密度算符随时间的变化。这种方法使用的波函数仍然是Kohn-Sham（KS）方程的波函数。该过程的核心在于无规相近似（random phase approximation, RPA）。
②含时密度泛函理论（TDDFT）:含时的微扰理论，它可以更好的描述体系对外场的响应函数χ(ω,t )。
在这个篇章中，我们将介绍第1种方法，在后续两个篇章中我们将介绍TDDFT方法。GW方法过于昂贵，而且其形式已经高于DFT的框架，在未来将会有专门的篇章介绍GW以及GW+BSE。
其中l 代表能带的指标，Ω 为使用的晶胞的体积。此处为没有使用大家熟悉的n作为能带指标道个歉，因为n 在前后文使用太多次了。其密度算符ρ0满足
现在我们给该体系施加一个随时间变化的外势Vi，体系会产生一个响应势Vs以屏蔽外场的扰动。此时体系的哈密顿量和密度算符会变为
由Poisson方程可知V和ρ1存在函数关系，为了简化计算，可以作线性化（linearize）处理，即将二者的乘积视为高阶小量，此时有
求和之后又对r做积分，这个过程非常繁琐，因此引入RPA：在电子浓度足够大时（固体中很容易满足这个条件），电子的位置在空间的分布是杂乱无规的，在由r决定的相位因子求和
这也可以理解为电子浓度足够高的情况下，相干叠加的数量级要远大于无规相位的指数相加。刚才的表述是经典动力学中的表述，所以才有“电子浓度”这类表述，源于David在1951年关于相互作用的电子气与等离激元的工作（Phys.Rev.82,625）中。笔者认为这种解释比较直观，所以在此使用。感兴趣的读者可以在李正中先生的《固体理论》的4-4中看到二次量子化表述下的量子运动方程的RPA，即将“异类”之间彼此激发产生的非线性项的求和归零，因此RPA也可以视为一种线性化操作。在微扰理论中，RPA又被称为环形图（ring diagram）近似或圈图修正，在此不继续展开讨论，想继续进阶的读者可以看：
体系对V 的响应表现在极化过程，极化会带来电荷密度的改变（束缚电荷），根据介质中的Maxwell方程组，电荷分为束缚电荷ρb和自由电荷ρf
从另一个角度理解上述关系，电荷密度的变化完全来源于为了屏蔽外势Vi而产生的极化。因此，体系响应外场而产生的屏蔽势Vs和电荷密度的变化可以有一个属于自己的Poisson方程
为了求解介电函数的矩阵，仍然需要繁琐的迭代，人们通常只考虑长波极限下的情况。长波极限对应的频率较高，比较符合可见光波段的吸收过程。当q远小于K时，可以近似的认为标量场f 在一个晶胞范围内的平均值就是标量场在K=0时的值，且f 的平均值的梯度就是f 梯度的平均值，于是
在只有垂直跃迁的情况下，l’ 即为导带CB，l 即为价带VB，再引入原子单位制，则虚部可以改写为
介电函数矩阵经常会存在不为0的非对角元，这是由于微观上单个偶极子处的电场强度并不是平均的宏观电场的强度，而应该是在单个电偶极子中心处的平均电场强度Eeff，这是一种局域场强。Eeff和宏观电场的差别来源正是介质中电偶极子的相互作用。在晶体的计算中，这种局域场效应则来源于周期性的势场。有多种方法可以修正局域场效应，我们将在今后的极化性质篇章，结合离子部分的介电函数等因素，详细讨论这个问题。
PWmat具有两个module可以利用RPA计算光学性质：module 18的前一部分和module 38。它们的前置工作是类似：
可以选择在这一步就将k-mesh取得很密，同时保证足够的空带数目。也可以选择在这一步使用较少的k-mesh点和空带数，然后在自洽计算的基础上做一个加密K点和空带数目的非自洽计算并输出波函数。
两个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 介绍
