关于HSE06杂化泛函计算能带的经验交流贴

duchenggong

各位老师、前辈、同学，最近在做HSE06的能带计算时查阅了很多资料和网站，但是总的来说没有一个系统的解析和教程，有些INCAR里面的参数设置也很不清楚，想在此开一个关于这方面的经验交流贴，也可以说是请教的帖子。

duchenggong

我首先说一下我做HSE06能带的步骤。刚接触理论计算不到一年，这也是自己查了一些资料总结出来的，不知道是否合理，希望各位前辈指点一下。我做的是2 x 2 x1 g-C3N4参照了部分vaspkit中文手册（比较老的一个版本的，应该是0.72的）
1、结构优化，

INCAR文件如下：
#### initial I/O ####
SYSTEM = 2x2-hg-C3N4 jiegou
ISTART = 1
ICHARG = 1
LWAVE = .FALSE.
LCHARG = .TRUE.
LVTOT = .FALSE.
LVHAR = .FALSE.
LELF = .FALSE.
# LORBIT = 11
# NEDOS = 2000

#### Ele Relax ####
ENCUT = 500 # 截断能，越高表示精度越高，默认400eV
ISMEAR = 0 # 和SIFGMA一起决定在做布里渊区积分时，如何计算分布函数。绝缘体0-0.05，金属1/2-0.2，态密度-5，SIGMA不设置
SIGMA = 0.05

EDIFF = 1E-5  # 控制自洽优化收敛的能量标准，即前后两次能量差小于这个值则收敛完成。默认1E-4
NELMIN = 5 #电子自洽的最小步数，默认为2
NELM = 300 #电子自洽的最大步数，默认60
GGA = PE #交换关联势的类型，VASP会根据POTCAR自动确定(我们用的赝势库所对应的为PE)，也可不用设置
LREAL = .FALSE.# 表面计算用Auto。决定投影算符在实空间还是倒空间计算，高精度计算用FALSE，在超过20个原子的系统中可以设置成Auto。
ISYM = 0 # 确定是否打开对称性，1、2、3为打开对称性，-1、0关闭对称性。
PREC = Normal # 控制计算精度的参数，会影响多个参数比如ENCUT。一般默认为Normal，高精度的计算比如声子计算时需要设置成Accurate。
IVDW = 10 # DFT-D2矫正范德华力。11为D3矫正。
# LDIPOL = .TRUE.# 两个控制表面偶极校正，对于上下表面不对称，会引入偶极，但是SCF收敛变慢
# IDIPOL = 3

#### Geo opt ####
EDIFFG = -0.001 # 设定结构优化的精度，正值表示能量，负值表示力，建议设置为负值（特别是声子计算时需要设置更小的值）
IBRION = 2 # 决定粒子如何运动。-1表示离子不动；0表示第一性原理分子动力学模拟；1表示准牛顿法进行结构优化；2表示共轭梯度法进行结构优化
POTIM = 0.2 # IBRION=0时，为时间步长，单位为fs；在结构优化计算时(1\2\3)，为力的缩放因子，默认为0.5
NSW = 300 # 离子运动的最大步数，在做结构优化和分子动力学模拟中设置。默认值为0，即离子不动。
ISIF = 2 # 0-不计算应力张量；1-结构优化时只计算总的应力，只优化离子位置。2~7计算应力张量。2只让离子位置做优化，晶胞矢量不做优化；3：同时优化离子位置、原胞形状和体积
ALGO=VeryFast #电子最优化的算法。

#### Mag ####
# ISPIN = 2 # 电子自旋极化的开关。默认为1，非磁性计算，自旋非极化体系；2.自旋极化体系。简单磁性体系直接设置成2，magmom不必设置。
# MAGMOM = 1*1.0 # 做磁性计算时的初始磁矩



2、PBE自洽，参照vaspkit手册，用vaspkit生成的INCAR，101，ST。这里我用的KPOINTS是Gamma Centered，5x5x1

INCAR文件如下：
Global Parameters
ISTART =  1            (Read existing wavefunction; if there)
# ISPIN =  2           (Spin polarised DFT)
# ICHARG =  11         (Non-self-consistent: GGA/LDA band structures)
LREAL  = .FALSE.       (Projection operators: automatic)
# ENCUT  =  400        (Cut-off energy for plane wave basis set, in eV)
PREC   =  Normal       (Precision level)
LWAVE  = .TRUE.        (Write WAVECAR or not)
LCHARG = .TRUE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT =             (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE.      (Real space distribution; supercells)
# NPAR   = 4           (Max is no. nodes; don't set for hybrids)
# NWRITE = 2           (Medium-level output)
# KPAR   = 2           (Divides k-grid into separate groups)
# NGX    = 500         (FFT grid mesh density for nice charge/potential plots)
# NGY    = 500         (FFT grid mesh density for nice charge/potential plots)
# NGZ    = 500         (FFT grid mesh density for nice charge/potential plots)

Static Calculation
ISMEAR =  0            (gaussian smearing method)
SIGMA  =  0.05         (please check the width of the smearing)
LORBIT =  11           (PAW radii for projected DOS)
NEDOS  =  2001         (DOSCAR points)
NELM   =  60           (Max electronic SCF steps)
EDIFF  =  1E-08        (SCF energy convergence; in eV)


这一步的目的是生成波函数和电子密度

3、HSE06自洽。正常应该是直接用杂化泛函做自洽，获取波函数和电荷密度，但是用PBE先得到电荷密度和波函数，再读取进行HSE06杂化泛函的自洽，速度会快一些。KPOINTS文件是Gamma
的3x3x1。INCAR是用vaspkit生成的，101，STH6（如果是D2矫正就在后面加上-D2），然后生成的INCAR文件如下：

Global Parameters
ISTART =  1            (Read existing wavefunction; if there)
# ISPIN =  2           (Spin polarised DFT)
# ICHARG =  11         (Non-self-consistent: GGA/LDA band structures)
LREAL  = .FALSE.       (Projection operators: automatic)
# ENCUT  =  400        (Cut-off energy for plane wave basis set, in eV)
PREC   =  Normal       (Precision level)
LWAVE  = .TRUE.        (Write WAVECAR or not)
LCHARG = .TRUE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT =             (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE.      (Real space distribution; supercells)
# NPAR   = 4           (Max is no. nodes; don't set for hybrids)
# NWRITE = 2           (Medium-level output)
# KPAR   = 2           (Divides k-grid into separate groups)
# NGX    = 500         (FFT grid mesh density for nice charge/potential plots)
# NGY    = 500         (FFT grid mesh density for nice charge/potential plots)
# NGZ    = 500         (FFT grid mesh density for nice charge/potential plots)

Static Calculation
ISMEAR =  0            (gaussian smearing method)
SIGMA  =  0.05         (please check the width of the smearing)
LORBIT =  11           (PAW radii for projected DOS)
NEDOS  =  2001         (DOSCAR points)
NELM   =  60           (Max electronic SCF steps)
EDIFF  =  1E-08        (SCF energy convergence; in eV)

HSE06 Calculation
LHFCALC= .TRUE.       (Activate HF)
AEXX   =  0.25        (25% HF exact exchange, adjusted this value to reproduce experimental band gap)
HFSCREEN= 0.2         (Switch to screened exchange; e.g. HSE06)
ALGO   =  ALL         (Electronic Minimisation Algorithm; ALGO=58)
TIME   =  0.4         (Timestep for IALGO5X)
PRECFOCK= N           (HF FFT grid)
! NKRED    = 2        (Reduce k-grid-even only, see also NKREDX, NKREDY and NKREDZ)
# HFLMAX   = 4        (HF cut-off: 4d, 6f)
# LDIAG    = .TRUE.   (Diagnolise Eigenvalues)


4、能带计算。K点文件用vaspkit生成，先用vaspkit的3，303，得到KPATH.in和PRIMCELL和HIGH_SYMMERY_POINTS。重启vaspkit，25，251，0.04，0.04，此时会根据KPATH.in 来生成KPOINTS文件。
我的INCAR文件如下：
#### initial I/O ####
SYSTEM = g-C3N4
ISTART = 1
ICHARG = 1
LWAVE = .FALSE.
LCHARG = .TRUE.
LVTOT = .FALSE.
LVHAR = .FALSE.
LELF = .FALSE.
LORBIT = 11
NEDOS = 2000

#### HSE06 ####
ALGO = DAMPED
TIME = 0.4
LMAXFOCK = 4
LHFCALC = .TRUE.
HFSCREEN = 0.2
PRECFOCK = FAST

#### Ele Relax ####
ENCUT = 500
ISMEAR = 0
SIGMA = 0.05
EDIFF = 1E-5
NELMIN = 5
NELM = 300
GGA = PE
LREAL = .FALSE.# 表面计算用Auto
ISYM = 0 # 对称性计算用0，
PREC = Nomarl
IVDW = 10
# LDIPOL = .TRUE.# 两个控制表面偶极校正，对于上下表面不对称，会引入偶极，但是SCF收敛变慢
# IDIPOL = 3

#### Geo opt ####
EDIFFG = -0.001
IBRION = -1
POTIM = 0.2
NSW = 0
ISIF = 2 # 2：只让原子位置做优化，晶胞矢量不做优化；3：晶胞优化，让晶格矢量优化

#### Mag ####
# ISPIN = 2 # 1.自旋非极化体系；2.自旋极化体系。简单磁性体系直接设置成2，magmom不必设置。
# MAGMOM = 1*1.0




问题1：第一步结构优化是否需要尽可能大的K点和尽可能高的截断能以及EDIFFG？尽可能大的话是多大？（我目前做的原子数不算多，60个原子）
问题2：我在计算能带的过程中的步骤是否合理？以及其中一些参数的设置？问题3：vaspkit生成的INCAR，好多参数（尤其是ICHARG）都没有读取，是用的默认值还是根据其他参数会自己调整？如果是它自己调整，那在不同情况下他自己调整完的值又是什么？
问题4：对于最后一步能带计算，此时的ICHARG是用 1 还是 11 ？杂化泛函的能带只能通过自洽来得到，不能和PBE的能带一样通过非自洽方法得到。但是我测试了一遍将ICHARG改成 1 的情况，速度上快了很多，而且数值上几乎一样，那用 1 而不用 11 的意义在哪？


以上只是我的一些拙见，欢迎各位前辈批评指正，也欢迎各位前辈在此帖讨论。

pwzhou

duchenggong 发表于 2022-10-19 19:19
我首先说一下我做HSE06能带的步骤。刚接触理论计算不到一年，这也是自己查了一些资料总结出来的，不知道是 ...

建议直接看vaspwiki上的说明：https://www.vasp.at/wiki/index.p ... _hybrid_functionals
HSE计算能带不需要读取电荷密度，可以提供PBE计算的波函数加速收敛，不提供也无所谓。

duchenggong

pwzhou 发表于 2022-10-19 21:08
建议直接看vaspwiki上的说明：https://www.vasp.at/wiki/index.php/Band-structure_calculation_using_hy ...

感谢老师，我先看一下wiki

duchenggong

请教一下对于建模时单层的来说，真空层的厚度会有什么影响？

我自己做了个测试，发现对band gap基本没有影响，但是会对价带顶和导带底以及费米能级的值产生线性影响

随着真空层的增大10→15，
VBM从-2.87变成了-3.77，
CBM从0.125变成了-0.76。
费米能级从-2.65变成-3.55。

请问这是什么原因？

Weldingspock

这么好的帖子竟然没人回，我留个名，现在也在算。之后分享一下经验

duchenggong

Weldingspock 发表于 2023-3-31 15:58
这么好的帖子竟然没人回，我留个名，现在也在算。之后分享一下经验

请问您最近进行的如何，我遇到一个提示是否可以讨论一下？
Your generating k-point grid is not commensurate to the symmetryof the lattice.This can causeslow convergence with respectto k-points for HF type calculations。

这个情况您遇见过吗？

Weldingspock

duchenggong 发表于 2023-4-12 10:24
请问您最近进行的如何，我遇到一个提示是否可以讨论一下？
Your generating k-point grid is not commen ...

没有遇到过，你的k点是自己写的吗。我都是用vaspkit生成的

Weldingspock

最近在老师的指导下算了一个HSE06能带的合作项目，在这里分享一下一些步骤。也希望大家帮忙指点改进，体系是CdTe cubic phase.
1. 首先是结构优化：
因为体系比较简单，我用了两种优化方式；
（1） ISIF=7直接优化，对于简单的体系可以这么用，INCAR如下：
Startparameters for structure optimization
# Gloabal paramters
ISTART = 0              # Read existing wavefunction =1
ICHARG = 2             # Default=2 if ISTART=0. Non-self-consistent: GGA/LDA band structures, =11
LWAVE =.FALSE.         # Write WAVECAR or not

# Parallelisation
LPLANE = .FALSE.       # Real space distribution; supercells
NCORE = 5              # Max is no. nodes; don't set for hybrids##
# NSIM = 1             # Divides k-grid into separate groups
KPAR = 2               # Divides k-grid into separate groups
# NPAR = 4             # Divides k-grid into separate groups

# Electronic relaxation
PREC = Normal          # Precision level, Normal or High
LREAL = .FALSE.           # Auto for molecule containing for than 20 ionics
ENCUT = 500            # Cut-off energy for plane wave basis set, in eV, default=max ENMAX in POTCAR,for bulk structure relaxation, ENCUT=1.3 * default
NELM = 60              # Max SCF steps
NELMIN = 5             # MIN SCF steps
EDIFF = 1E-05          # SCF energy convergence, this number is quite precise
ISPIN = 1              # Closed shell DFT

# Ionic relaxation
EDIFFG =  -0.01        # Ionic convergence; eV/AA^3
NSW = 1000             # MAX Inonic relaxation steps
ISIF = 7               # For single point enegy
IBRION = 2            # = -1 for NSW = -1 or 0
ISMEAR = 0             # Gaussian smearing; metals:1
SIGMA = 0.05           # Smearing value in eV; metals:0.2
（2）使用BM方程，改变POSCAR中的晶胞缩放系数，计算单点能。可以参考Ex33 晶格参数的确定（Birch-Murnaghan状态方程-1） | Learn VASP The Hard Way (bigbrosci.com)。 将计算出来的单点能进行拟合（二次曲线），找到能量最低点。在本案例中，这种方法比（1）中计算出来的晶格参数更加稳定，能量更小。但是也需要注意，PBE泛函方法优化的晶格参数往往会过大或者过小的偏离实验值，这个是正常的。这里贴一个单点能计算的INCAR，晶格缩放系数为1.0
Startparameters for structure optimization
# Gloabal paramters
ISTART = 0              # Read existing wavefunction =1
ICHARG = 2             # Default=2 if ISTART=0. Non-self-consistent: GGA/LDA band structures, =11
LWAVE =.FALSE.         # Write WAVECAR or not

# Parallelisation
LPLANE = .FALSE.       # Real space distribution; supercells
NCORE = 5              # Max is no. nodes; don't set for hybrids##
# NSIM = 1             # Divides k-grid into separate groups
KPAR = 2               # Divides k-grid into separate groups
# NPAR = 4             # Divides k-grid into separate groups

# Electronic relaxation
PREC = Normal          # Precision level, Normal or High
LREAL = .FALSE.           # Auto for molecule containing for than 20 ionics
ENCUT = 500            # Cut-off energy for plane wave basis set, in eV, default=max ENMAX in POTCAR,for bulk structure relaxation, ENCUT=1.3 * default
NELM = 60              # Max SCF steps
NELMIN = 5             # MIN SCF steps
EDIFF = 1E-05          # SCF energy convergence, this number is quite precise
ISPIN = 1              # Closed shell DFT

# Ionic relaxation
EDIFFG =  -0.01        # Ionic convergence; eV/AA^3
NSW = 0                   # MAX Inonic relaxation steps
ISIF = 2               # For single point enegy
IBRION = -1            # = -1 for NSW = -1 or 0
ISMEAR = 0             # Gaussian smearing; metals:1
SIGMA = 0.05           # Smearing value in eV; metals:0.2


2. PBE泛函进行自洽计算
得到WAVECAR，为之后的HSE能带计算做准备。（为什么不进行HSE自洽楼上已经说得很清楚了）

3. 确定HSE中的AEXX参数，即确定HF泛函在hybird function中所占的比例。
vasp中默认AEXX=0.25，但是这个数字计算出来的band gap往往和实验值对不上。所以就要像1中一样，调控AEXX的大小，使算出来的band gap和实验值一致。计算band gap最好使用原胞，原子少。一般band gap与AEXX的关系是线性的，大家只需要找到合适的值就好。在这里贴一个hse计算能带的INCAR，请注意要读取上一步的wavecar。
# Gloabal paramters
ISTART = 1             # Read existing wavefunction =1
ICHARG = 1             # Default=2 if ISTART=0. Non-self-consistent: GGA/LDA band structures, =11
LWAVE =.FALSE.         # Write WAVECAR or not

# Parallelisation
LPLANE = .FALSE.       # Real space distribution; supercells
NCORE = 4              # Max is no. nodes; don't set for hybrids##
# NSIM = 1             # Divides k-grid into separate groups
KPAR = 4               # Divides k-grid into separate groups
# NPAR = 4             # Divides k-grid into separate groups

# Electronic relaxation
PREC = Normal          # Precision level, Normal or High
# LREAL = Auto         # Auto for molecule containing for than 20 ionics
ENCUT = 500            # Cut-off energy for plane wave basis set, in eV, default=max ENMAX in POTCAR,for bulk structure relaxation, ENCUT=1.3 * default
NELM = 60              # Max SCF steps
NELMIN = 5             # MIN SCF steps
EDIFF = 1E-05          # SCF energy convergence, this number is quite precise
ISPIN = 1              # Closed shell DFT

# Ionic relaxation
# EDIFFG =  -0.02        # Ionic convergence; eV/AA^3
NSW = 0                # MAX Inonic relaxation steps
ISIF = 2               # For single point enegy
IBRION = -1             # = -1 for NSW = -1 or 0
ISMEAR = 0             # Gaussian smearing; metals:1
SIGMA = 0.05           # Smearing value in eV; metals:0.2

# Hybrid Functional
LHFCALC = .TRUE.
HFSCREEN = 0.2
AEXX = 0.25
TIME = 0.4
PRECFOCK = Fast


4. 经过3中的筛选得到合适的AEXX，接下来对于这个体系的一些其他变化（例如晶胞膨胀缩小对band gap的影响），也可以使用同样的AEXX。

duchenggong

Weldingspock 发表于 2023-4-12 15:51
没有遇到过，你的k点是自己写的吗。我都是用vaspkit生成的

就是vaspkit 302，然后vaspkit 251，用的Gamma，0.04 0.04。都出现这个提示了，但是基本上都可以正常收敛，就是有一个结构200多步都不收敛，用了35＋个小时。

duchenggong

Weldingspock 发表于 2023-4-12 19:25
最近在老师的指导下算了一个HSE06能带的合作项目，在这里分享一下一些步骤。也希望大家帮忙指点改进，体系 ...

请问体系大概多少原子呢？

Weldingspock

duchenggong 发表于 2023-4-12 20:38
请问体系大概多少原子呢？

我这个用原胞算的，只有两个原子

Weldingspock

duchenggong 发表于 2023-4-12 20:35
就是vaspkit 302，然后vaspkit 251，用的Gamma，0.04 0.04。都出现这个提示了，但是基本上都可以正常收敛 ...

有可能是你用vaspkit生成的高对称点有问题（K-PATH in），建议用material studio的布里渊区功能找

duchenggong

Weldingspock 发表于 2023-4-13 12:22
有可能是你用vaspkit生成的高对称点有问题（K-PATH in），建议用material studio的布里渊区功能找

好的，谢谢。还有个问题想请教一下。我看您INCAR中的ISMEAR=0，计算完能带之后，可以用vaspkit直接输出态密度，请问这个态密度可以直接用吗？或者说这个不能用，得需要重新提交专门计算DOS的任务呢？

Weldingspock

duchenggong 发表于 2023-5-15 14:07
好的，谢谢。还有个问题想请教一下。我看您INCAR中的ISMEAR=0，计算完能带之后，可以用vaspkit直接输出态 ...

不能用，要专门提交计算DOS的任务。我记得dos的话ISMEAR应该设置成-5