TDDFT计算软件Octopus学习笔记(四)：带隙优化方法

密度泛函理论是一种基态的理论，然而能带结构和带隙属于物质的激发态特性。因此DFT的固有缺点就是会低估带隙能量。这一缺点可以通过不同的方式加以修正。本文介绍两种方法。

(一) DFT+U方法 (以ZnO为例)

氧化锌ZnO的实验带隙能量为3.3eV。 若如学习笔记(三)中，设：

PseudopotentialSet = pseudodojo_pbe

则Octopus 会使用GGA泛函，所得能带结构如下：

带隙能量只有0.79eV，显著低于实验值3.3eV。因此需要进行修正。

DFT+U的基本思路是将局域电子之间的强库伦相互作用写为一项类Hubbard项。 例如，对于ZnO这种过度金属氧化物来说，Zn的电子排布是 $1s^{2}2s^{2}2p^{6}3s^{2}s{p}^{6}3d^{1}04s^2$, 局域电子例如3d电子轨道之间彼此会有很大的重叠，3d电子将会彼此靠近，这将会增强库伦相互作用项到不容忽视的地步，电子关联变得十分重要。而LDA和GGA无法正确描述这种情况。
DFT+U模型将这种相互作用以参数U和J进行描述，这种参数通常是半经验的，这里我们以文献 Phys. Rev. X 5, 011006(2015) 为参考，选定 对于Zn的3d轨道U=12.8eV，O的2p轨道U=5.29eV。
最终，DFT+U修正后的ZnO带隙能量为2.79eV， 虽然仍小于实验值3.3eV，但是相对于GGA已经有了很大改进。

DFTULevel = dft_u_empirical               # 使用DFT+u 来扩展带隙
%Species
 "Zn" | species_pseudo | hubbard_l | 2 | hubbard_u | 12.8*eV
 "O"  | species_pseudo | hubbard_l | 1 | hubbard_u | 5.29*eV
%

代码注释： DFTULevel = dft_u_empirical 指定使用经验参数的DFT+U模型进行计算，
而 DFTULevel =dft_u_none 可设定不使用DFT+U 方法。
hubbard_l 指明 +U 所应用的轨道， 0,1,2,3 即 s,p,d,f 轨道。
hubbard_u 表明所 +U 能量的大小。

(二) metaGGA泛函方法(以金刚石©为例)

同样，金刚石的实验带隙大小约为5.4eV，如果使用 LDA泛函，则 得到带隙仅为 4.2eV， 因此，使用更高阶的泛函 metaGGA 加以修正。 我们将要使用的metaGGA 叫 TB09， 是一种轨道依赖的泛函，它依赖于动能密度，然而在重启信息文件夹 restart 中并没有写入相应信息，所以 按照 上篇 所言的方法会报出致命错误(fatal error).

**************************** FATAL ERROR *****************************
*** Fatal Error (description follows)
*--------------------------------------------------------------------
* For an orbital-dependent functional, all occupied orbitals must be provided.
* Not all the occupied orbitals could be read.
* Please run a ground-state calculation first!
**********************************************************************

这意味这如果我们想要读取gs计算的重启信息来运行unocc计算，这两者需要有相同的波函数， 即k点的数量和态相同。
之前我们计算 gs 基态时，所设k为 ：

%KPointsGrid
  12  |  12  |  7
%
KPointsUseSymmetries = yes

然而进行非占据态unocc计算时，所设k为：

%KPointsPath        #k点扫描路径,见“注1”
  14 |  10 | 20 | 14 | 14 | 10    #标明各个点之间的路径需要取样几个k点
 0.0 | 0.0 | 0.0  # Gamma point  #所以有7个点的坐标，只有6段路径需要标明取样点数
 0.0 | 0.5 | 0.0  # M  point
 -1/3| 2/3 | 0.0  # K point
 0.0 | 0.0 | 0.0  # Gamma point
 0.0 | 0.0 | 0.5  # A point
 0.0 | 0.5 | 0.5  # L point
 -1/3| 2/3 | 0.5  # H point
%

KPointsUseSymmetries = no    #停止使用对称性

这就是报错的根源。
因此，对于gs计算，应同时加入两种k点取样方法，inp文件应该改为：

CalculationMode =  gs
ExperimentalFeatures = yes
FromScratch = yes
PeriodicDimensions = 3
BoxShape = parallelepiped
ParKPoints = auto
XCFunctional = mgga_x_tb09 + lda_c_pw    # 交换部分设为TB09的metaGGA， 关联部分仍用LDA
FilterPotentials = filter_none                   # 使用势屏蔽，若不屏蔽，该泛函对应的 gs 计算很难收敛

a = 3.57*angstrom  #2.527/sqrt(2)
Spacing = 0.2

%LatticeParameters
  a | a | a
%

%LatticeVectors
 0.0 | 0.5 | 0.5
 0.5 | 0.0 | 0.5
 0.5 | 0.5 | 0.0
%
%ReducedCoordinates
 "C" |   0.0    | 0.0     | 0.0
 "C" |   3/4    | 3/4     | 3/4
%
nk = 16
%KpointsGrid
 nk | nk | nk
%
%KPointsPath
  20 |  20 |  30
 0.5 | 0.0 | 0.0  # L point
 0.0 | 0.0 | 0.0  # Gamma point
 0.0 | 0.5 | 0.5  # X point
 1.0 | 1.0 | 1.0  # Another Gamma point
%

KPointsUseSymmetries = no

这样， unocc计算的inp文件为：

CalculationMode =  unocc
ExperimentalFeatures = yes
FromScratch = yes
PeriodicDimensions = 3
BoxShape = parallelepiped
ParKPoints = auto
XCFunctional = mgga_x_tb09 + lda_c_pw
FilterPotentials = filter_none

a = 3.57*angstrom  #2.527/sqrt(2)
Spacing = 0.2

%LatticeParameters
  a | a | a
%

%LatticeVectors
 0.0 | 0.5 | 0.5
 0.5 | 0.0 | 0.5
 0.5 | 0.5 | 0.0
%
%ReducedCoordinates
 "C" |   0.0    | 0.0     | 0.0
 "C" |   3/4    | 3/4     | 3/4
%
nk = 16
%KpointsGrid
 nk | nk | nk
%

%KPointsPath
  20 |  20 |  30
 0.5 | 0.0 | 0.0  # L point
 0.0 | 0.0 | 0.0  # Gamma point
 0.0 | 0.5 | 0.5  # X point
 1.0 | 1.0 | 1.0  # Another Gamma point
%
KPointsUseSymmetries = no
ExtraStates = 10
ExtraStatesToConverge = 8
output = dos

最终所得带隙能量为4.96eV，比LDA效果好很多。
该方法不只对metaGGA生效，其他的杂化泛函，如果是轨道依赖的泛函，都需要注意这一点。