3. 光学性质的计算

• 3.1 用epsilon.x计算光学性质

• 3.2 用TDDFPT模块计算光学性质

• 3.3 静态及光学介电常数计算

3.1 用epsilon.x计算光学性质

这里，首先介绍一种基于LDA/GGA的通过Kohn-Sham轨道跃迁矩阵元计算介电常数的QE后处理程序epsilon.x。编译了QE的PP模块（make pp）之后，目录qe/PP/src下会生成这个可执行文件。epsilon.x计算能级间偶极跃迁得到介电常数，根据以下（固体中的费米黄金定则）^[1]得到介电常数虚部[7]：

而介电常数实部可以由虚部通过Kramers-Kronig关系变换得到。通过能带结构，可以计算得到全部的光学常数。

以计算闪锌矿结构立方氮化硼的介电函数为例。epsilon.x计算目前仅支持模守恒赝势。首先优化晶格常数，运行pw.x < vc.inp > vc.out，vc.inp如下：

&CONTROL calculation='vc-relax', disk_io='low', prefix='pwscf', pseudo_dir='./', outdir='./tmp', verbosity='high' tprnfor=.true., tstress=.true., forc_conv_thr=1.0d-5/&SYSTEM ibrav= 0, nat= 2, ntyp= 2, occupations = 'fixed' ecutwfc= 80, ecutrho = 320,/&ELECTRONS electron_maxstep = 100 conv_thr = 1.0d-9 mixing_mode = 'plain' mixing_beta = 0.8d0 diagonalization = 'david'/&IONS ion_dynamics='bfgs'/&CELL press_conv_thr=0.1/ATOMIC_SPECIES B 10.811 B_ONCV_PBE-1.0.upf  N 14.00674 N.oncvpsp.upfCELL_PARAMETERS (angstrom) 1.807500000 1.807500000 -0.000000000 0.000000000 1.807500000 1.807500000 1.807500000 -0.000000000 1.807500000ATOMIC_POSITIONS (crystal)B 0.000000000 0.000000000 0.000000000N 0.250000000 0.250000000 0.250000000K_POINTS {automatic} 13 13 13 0 0 0

赝势文件见 B_ONCV_PBE-1.0.upf^[2]N.oncvpsp.upf^[3]。

得到晶格常数PBE计算值3.6215?。光学性质需要较密的k点，并且需要关掉对称性，设置nosym=.true,noinv=.true，计算的带个数nbnd要大一些，做一次nscf计算：

&CONTROL calculation='nscf', disk_io='low', prefix='pwscf', pseudo_dir='./', outdir='./tmp', verbosity='high' tprnfor=.true., tstress=.true., forc_conv_thr=1.0d-5/&SYSTEM ibrav= 0, nat= 2, ntyp= 2, occupations = 'smearing', smearing='gauss', degauss=1d-9, ecutwfc= 80, ecutrho = 320, nbnd=30 nosym=.true, noinv=.true,/&ELECTRONS electron_maxstep = 100 conv_thr = 1.0d-9 mixing_mode = 'plain' mixing_beta = 0.8d0 diagonalization = 'david'/&IONS ion_dynamics='bfgs'/&CELL press_conv_thr=0.1/ATOMIC_SPECIES B 10.811 B_ONCV_PBE-1.0.upf  N 14.00674 N.oncvpsp.upfCELL_PARAMETERS (angstrom) 1.810733563 1.810733563 0.000000000 -0.000000000 1.810733563 1.810733563 1.810733563 -0.000000000 1.810733563ATOMIC_POSITIONS (crystal)B 0.000000000 0.000000000 0.000000000N 0.250000000 0.250000000 0.250000000K_POINTS {automatic} 13 13 13 0 0 0

以下内容保存为epsilon.inp，运行epsilon.x < epsilon.inp > epsilon.out：

&inputpp outdir='./tmp' calculation='eps'/&energy_grid smeartype='gauss' intersmear=0.50 intrasmear=0.0 wmin=0.0 wmax=60.0 nbndmin=1 nbndmax=0 nw=2000 shift=0.0/

输出4个文件eels_pwscf.dat,epsi_pwscf.dat,epsr_pwscf.dat,ieps_pwscf.dat，分别是eels谱（电子能量损失谱），介电函数的虚部、实部，虚轴上的介电函数。前三个结果见下图，与实验及其他计算结果比较见[8]。

固体的线性光学性质，介电常数、电导率、折射率一般来说都是复数，知道其一可以确定其余两个，进而还可以得到吸收系数，反射率，三者关系如下表[9]：

吸收系数

反射率

3.2 用TDDFPT模块计算光学性质

计算固体光学性质的方法，比较先进的是通过GW-BSE方程加入电子的自能（GW）、电子-空穴相互作用修正（BSE），以及含时密度泛函理论(TDDFT)等方法[17]，在QE中有TDDFPT，GWL，Yambo等相关模块。这里介绍一种计算量相对比较小的TDDFPT方法。

TDDFT与GW-BSE计算相比，GW-BSE计算使用格林函数，直接构造了物理过程，是计算光学性质的标准方法；TDDFT基于电子密度，包含了各种作用但并不那么直观，重要的是TDDFT计算量比GW-BSE小很多。这两种方法是互补的关系，原则上都是精确的。TDDFT使用的交换关联泛函十分关键，原则上需要非局域交换关联近似，如Hybrid functional、meta GGA、OEP等（后两种暂未实现，hybrid只支持turbo_lanczos.x）。

TDDFPT模块[15]，包括适用于分子体系的turbo_lanczos.x（以及turbo_davidson.x），和计算固体的turbo_eels.x，以及它们的后处理程序turbo_spectrum.x。turbo_lanczos.x只支持单个k点（ 点），对于固体要用超胞近似。下面使用turbo_eels.x模块，turbo_eels.x需要指定转移动量，一般地，介电函数是动量和频率的函数，这在电子能量损失谱（EELS）中较容易理解，介电函数动量谱的意义是体系在空间非均匀电场下的响应，EELS实验上是入射、出射电子的动量变化，光学性质是动量趋向于0的极限。

同3.1节的BN材料，首先进行一次scf计算。

&CONTROL calculation='scf', disk_io='low', prefix='pwscf', pseudo_dir='./', outdir='./tmp', verbosity='high' tprnfor=.true., tstress=.true., forc_conv_thr=1.0d-5/&SYSTEM ibrav= 0, celldm(1) = 1.8897261328856432, ! a.u. to Angst nat= 2, ntyp= 2, occupations = 'smearing', smearing='gauss', degauss=1d-9, !occupations = 'fixed', ecutwfc= 80, ecutrho = 320, nbnd=30,/&ELECTRONS electron_maxstep = 100 conv_thr = 1.0d-9 mixing_mode = 'plain' mixing_beta = 0.8d0 diagonalization = 'david'/&IONS ion_dynamics='bfgs'/&CELL press_conv_thr=0.1/ATOMIC_SPECIES B 10.811 B_ONCV_PBE-1.0.upf N 14.00674 N.oncvpsp.upfCELL_PARAMETERS (alat= 1.88972613) 1.810733563 1.810733563 0.000000000 -0.000000000 1.810733563 1.810733563 1.810733563 -0.000000000 1.810733563ATOMIC_POSITIONS (crystal)B 0.000000000 0.000000000 0.000000000N 0.250000000 0.250000000 0.250000000K_POINTS {automatic} 8 8 8 0 0 0

turbo_eels.x的输入如下，需要指定一个接近0的动量q值[16]。itermax是时间步数。

&lr_input prefix='pwscf', outdir='./tmp' restart_step=10, restart=.false./&lr_control itermax = 5000, ipol=1, q1 = 0.001d0, q2 = 0.0d0, q3 = 0.0d0,/

运行turbo_spectrum.x，输入如下，注意使用了插值extrapolation，展宽值epsil选取对结果影响较大。

&lr_input prefix='pwscf', outdir='./tmp' eels = .true. itermax0 = 5000 itermax = 10000 extrapolation="osc" epsil=0.0350 units=1 start=0.0d0 increment=0.001d0 end=60d0/

输出pwscf.plot_eps.dat，画图如下。

3.3 静态及光学介电常数计算

3.3.1 加电场的自洽计算

对于slab模型、分子等，有时由于对称性的破坏，在平面波程序的周期性边界条件下，会造成人为的电场，需要加入一个dipole去消除；另外，有时需要研究体系在电场下的电荷变化，如3.3.2节中的介电常数计算，需要加电场的计算。在真空层中加入一个电偶极子，也就是引入一个沿着z方向（z是垂直表面方向，由edir设置）的锯齿形势能，这个势能在slab处上升，上升的梯度为电场强度eamp，在真空一个很小的宽度（eopreg）内下降。计算是否要做relax，根据实际情况定。设置需要用到edir,emaxpos,eopreg,eamp,tefield和dipfield。当slab位于Cell中间（约0.5）时（输入其余部分略）：

&CONTROL tefield=true., dipfield=.true./&SYSTEM edir = 3, eamp = 0.001, emaxpos=0.99, eopreg=0.01,/

偶极矩的位置在z方向的0.99~1.0之间。

QE中另一种加电场的方法是用berry phase（BP）加入电场，与dipole电场适用的场合不同，BP方法要求体系是没有真空层的，或者有真空层但是电场在slab的面内（x，y方向）。BP电场设置参数

 lelfield=.true., nberrycyc=3... efield_cart(3)=0.001d0

注意：对于导体不适用加电场，（1）dipole方法，金属slab加电场会让正负电荷在两个表面积累，这种电荷转移的量很大，电荷自洽计算几乎不可能收敛。（2）berry phase方法，由于berry phase是在倒空间一个闭合的曲面上的积分，金属的占据态对应k点集合不是闭合的曲面，无法定义Berry phase；对于绝缘体，加的电场过大会将体系带隙变小直至成为导体而无法计数。

3.3.2 静态及光学介电常数计算

在3.3.1节的基础上，可以计算固体的静态及光学介电常数，参考安装包中的q-e-qe-6.5\PW\examples\example10，原理和参考文献看一下其中的README，分别做两次相反方向小电场的自洽计算，得到极化的改变量，从而得到介电常数。

计算静态介电常数，只要在上面光学介电常数的基础上，允许结构做relax就可以得到。

在example10中，使用的是berry phase电场，对于含有真空的cell，则可以用dipole电场。两种加电场方法都可以得到比较准确的结果。需要注意的是相关物理量（总能、力）的收敛问题，相关阈值要适当调小，经过测试，即使非常接近铁电的体系仍然可以得到合理的结果。

参考文献：

[7] F. Bassani, G. Parravicini, Electronic states and optical transitions in solids. Pergamon Press (1975).

[8] Anna Tararan et al. Optical gap and optically active intragap defects in cubic BN. Phys. Rev. B 98, 094106 (2018).

[9] Martin Dressel, George Gruner, Electrodynamics of Solids-Optical Properties of Electrons in Matter. Cambridge University Press (2002).

[15] I. Timrov, N. Vast, R. Gebauer, S. Baroni, turboEELS — A code for the simulation of the electron energy loss and inelastic X-ray scattering spectra using the Liouville–Lanczos approach to time-dependent density-functional perturbation theory, Computer Physics Communications (2015).

[16] Iurii Timrov, Maxime Markov, Tommaso Gorni, Michèle Raynaud, Oleksandr Motornyi, Ralph Gebauer, Stefano Baroni, and Nathalie Vast. Ab initio study of electron energy loss spectra of bulk bismuth up to 100 eV. Phys. Rev. B 95, 094301 (2017).

[17] Onida, Giovanni and Reining, Lucia and Rubio, Angel. Electronic excitations: density-functional versus many-body Green’s-function approaches. Rev. Mod. Phys., 74, 601(2002).

参考资料

[1]

公式: https://en.wikipedia.org/wiki/Fermi%27s_golden_rule

[2]

B_ONCV_PBE-1.0.upf: https://yyyu200.github.io/DFTbook/img/B_ONCV_PBE-1.0.upf

[3]

N.oncvpsp.upf: https://yyyu200.github.io/DFTbook/img/N.oncvpsp.upf